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Abstract 

— 

_ A first order perturbation analysis is developed for the 


restricted n-body problem of interplanetary flight. The 
4 linearized perturbation equations, derived from a generalized 

functional representation of the first and second integrals of 
- the dynamical differential equations, appear in a modified 
* _ tensor notation (without recourse to the summation conven- 
tion) which underscores their inherent structural simplicity. 
At the heart of this formulation is a time dependent Green’s 
- function, or influence tensor, which is derived for both fixed 
_ terminal times (upper limits of integration) and open terminal 
: times defined by some terminal constraint. The text concludes 
_ with a more or less complete mathematical development of 
two techniques for evaluating the fixed terminal time Green’s 
function. Direct integration, which constitutes the first 
method, is at once the most straightforward and obvious of 
the pair. It is included here simply for contrast with the second 
method, an indirect form of “integration” that utilizes the 
_ properties of a system of equations adjoint to a linearized 
form of the trajectory differential equations. A self-contained 

mathematical derivation of the latter technique appears, for 
convenience, in an appended discussion. 


_ Introduction 


Of late, analyses designed to predict the optimum 
magnitude and spacing of corrective thrust application 
for maneuverable space vehicles have solicited con- 
siderable interest. At the same time, the companion 
problem of trajectory and orbit determination in the 
face of measurement errors inherent in available track- 
ing systems deserves, and has received, attention. A 
common requirement for both types of studies is the 
capability to predict a vehicle’s departure from a so- 

called precalculated nominal or standard trajectory 
at a specified point in time, caused by position and/or 
velocity increments introduced at a previous time. 
This capability is most generally achieved through a 
mathematical formulation based upon certain linearity 
assumptions which require only that the perturbations 
considered (be they actual velocity increments caused 
by thrust impulses, or position and/or velocity measure- 
ment errors) are small quantities. The linearized ap- 
proach ultimately involves the calculation of a Green’s 
function whose components are sometimes referred to 
as error coefficients. 


1 Manuscript received December 2, 1960, revised August 
4, 1961. 

2 Research Engineer, Grumman Aircraft Engineering 
Corporation, Bethpage, New York. 


Green’s Functions for Space Trajectory 
Perturbation Analysis 


Joseph C. Dunn? 


The primary objective of this paper is to put before 
the analyst a brief but unified theoretical presentation 
of linearized perturbation theory as it applies to inter- 
planetary trajectory studies. An attempt has been 
made to place the emphasis, as much as is possible, on 
mathematical rigor in order to complement the exposi- 
tions of practical application which have already been 
made in some detail by Battin, Breakwell, Laning, and 
Pfeiffer, {1], [2], [3], [4], and [5]. Where the need arises, 
a few simple examples have been inserted for clarifica- 
tion. 


Linearized Trajectory Perturbation Equations 
for the Restricted Three-Dimensional n- 


Body Problem 


The differential equations describing the motion 
of several mutually attractive central gravitational 
force fields which are in transit with respect to an 
inertial space,’ may be written as follows [6]: 


te = =I ym 7 see 2s) ) (1a) 
j=1 aj 
s = (yi; — y;) 4=1n 
= Yo Pease 1b 
us K Dim, ee * j=l n3j At (1b) 
Aver coh fs ed | (1e) 
j=l Pea 


riz = [Ca ") — xj)” + (yi — yi)” + (2: — 2;)'] he (1d) 


where 2; , y;, and 2; define the position of the 7th mass 
point relative to an inertial cartesian space, and the 
superscripted dot indicates a single differentiation in 
time. 

Equations (1a), (1b), and (le) may be expressed 
in a symbolic form which simply states the functional 
dependence, e.g.: 


&; = fit, Ys, 23), (2a) 

Yi = Gil%i, Ys , 23), t=1,n;j =1, (2b) 
and 

Be = hi 2; , Ys» 2;)- (2c) 


At time ¢ = 1, define a set of nominal initial conditions 
in a position-velocity configuration space such that: 


Qi. = yl); b=-1) 6397 ="; 


3 General n-body problem of celestial mechanics. 
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where gj = 235 932 = Us 903 = 27, 14 = Tj, Qs = Ys» 
gis = 2;- Then the corresponding first and second 
integrals of Eqs. (2a), (2b), (2c) are seen to depend 
solely upon these initial conditions and a nominal 
terminal time, ¢; , (upper limit of integration) counted 
from 7, 1.€:, 


qiz (ty) = Fi(qylr), aks 
= 1, 630 Le. 


(3) 


A question presents itself: What increments, Aqix , 
in the nominal variables qiz(t;) are produced by the 
perturbations Agqj;(r) and At;, introduced at time 
7 < t,? Or, from a physical point of view, how will a 
number of small variations in individual component 
orientations and velocities introduced at time 7, 
ultimately affect the subsequent n-body configuration 
at a new terminal time ¢; = t; + At;? The answer, of 
course, can be arrived at by simply reintegrating the 
n-body equations of motion for the new set of initial 
conditions gji(7) + Agji(7) and extending the upper 
limit of integration to ¢; + At,;. But, if the perturba- 
tions Aqji(7) and At; are constrained to remain ‘“‘small,”’ 
a more convenient approach to the problem is justified. 

By way of definition: 


Agqin(ts) = Gie(ts) — qie(ty), (4) 


where the g;,(ty) are nominal (unperturbed) param- 
eters evaluated at the nominal terminal time t; and 
the @i(¢;) are off-nominal parameters evaluated at 
new terminal time tr; = ty; + At; , following the in- 
troduction of Agqj:(7) perturbations at time 7. The 
Gix(ty) may be estimated, in first approximation, by 
means of a truncated Taylor series expansion of Kq. 
(3) in the neighborhood of q;,(7) and t;. Then, to the 
exclusion of second and higher order terms; 


n 6 
= z OF ix 
Gir(Ey) = gir(ty) + >) >> — (7) Aqi(7) 
jai lai O51 
= (5) 
ik 
SEE AT 
“I at ( 1) Hee) 
assuming that the derivatives a and ae exist and 
qjl ) 


are continuous in the neighborhood of the expansion. 
With the aid of Eqs. (4) and (5) one finds that: 


n 


6 
7 OF OF... 
Moms) 52 5s oan een ere 
qik (ly 2, rs, Mp T)Agy\7) +e a (ty YAR 6G) 
Also, consideration of Eq. (3) reveals that OF ix eee 
Ot 
d . . 
di qik(tz) = Gir(ty) since; 
ne Ps = (4) =n nee 4d # 
40K \ of — Yir¥ = =~ 5 
dt ‘ jah [a OQe7 See qit(r) (7) 
OF 5, 
i Ss: t 
zs dt (dy) 


i fa pus 4 ‘ 
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and; 
d . 
at 2" 

Equation (6) can now be written in a somewhat 
more descriptive notation, viz: | 


n 6 4 ne ; { 
Agi(ty) = Dae — (r)Ag(r) + Gin(ts)Aty ( | 


a t=2 


(the initial conditions are not (8 
functions of time). 5 


(rr 0; 


which comprises the set of 6n linearized perturbatio 4) 
equations in 6n unknowns, Agi(t;), for the general 
n-body problem. Now consider its special applications 
to a restricted n-body problem. 3} 

Suppose that only the mth particle of all the con+ 
stituent elements in an n-body group is subject tot 
perturbations, 1.e.; 1 


Agyi(r) = 8jm Agylt); Sim = ° 


Then Eq. (9) simplifies to: 


6 


s OF; 
Agu (ty) = pa, : 


i= Oat 


(7) Aqmilr) + Gir (ts)Aty - (10) 


If it is further assumed that the same mth particle? 
exerts an insignificant influence on the remaining 
n —-1 bodies, then it follows that: 


ee (7) = Og OP ix (7); Oim = ee m= m i 


Odmt Odmt1 i =m i 
Consequently ; 
> *. aF 
Adm (ty) = a 3 = (1) Agmi(t) + Gm (ty) Aty (11a) 
= ml 
and, 
Aqi(ts) = qx(ty) At, ; 1#m. (11b) 


It is significant that the assumptions which make 
Eqs. (11a) and (11b) a special case of Eq. (9) have 
real physical validity for interplanetary trajectory 
perturbation mechanics. Most certainly, a spacecraft 
exerts a totally negligible influence on planetary and/or 
stellar elements in its vicinity: also, one commonly 
assumes a perfect knowledge of planetary kinematics 
and dynamics when computing trajectories. Conse- 
quently, in this light, the spacecraft is the only member 
of the n body group (planets + Sun + spacecraft) 
which is not exempt from perturbative influence. 
It then must follow that changes in a spacecraft’s 
“state” (position and velocity) at a terminal point 
may only be attributed to deviations from its initial 
state at a previous time. At this point, Eqs. (lla) are. 
immediately applicable provided that the deviations 
from the standard initial state are “small.” 

In actuality, planetary kinematics and dynamics are 
not known precisely, and the uncertainties which do. 
exist (e.g., in the astronomical constants, ete.) can re 
sult ina definite dispersion from the precalculated termi- 
nal state even if initial conditions are met exactly. 


An assessment of the influence of uncertainties of this 

ort can be made by elevating the parameters in doubt 
to the level of a state variable, and then proceeding to 
construct a set of expanded perturbation equations. 
‘The procedure is identical to that described in earlier 
pages and, for example, might start from a set of dif- 
ferential equations of the following sort: 


é; = fia; , Kr) 
K, — 0 


z where the K;’s are here taken to be coefficients 
the dynamical equations whose values are in doubt. 
Unfortunately, further exploration of this trend of 

‘thought, however fruitful, does not lie within the pro- 
scope of this paper: subsequent developments 

ill therefore begin with Eq. (lla) and the implica- 
tions contained therein. 

Henceforth, the subscript m in Eq. (11a) will be 
Biicleted and it is to be understood that the parameters, 
q, are to be associated only with vehicle coordinates 
and coordinate velocities. Then Eq. (1la) goes over 
to: 


: 
ee 


a teh. Th, 


6 
alt) = oS (rags) + de(t)ate* (12) 
where, in summary: Aq,(t;) is a vector whose com- 
ponents are incremental displacement and _ velocity 
deviations, at time 7; = ¢t; + At; , from nominal trajec- 
tory variables at time ty; Aqgi(7) is a vector whose 
components are displacement and velocity perturba- 
tions introduced at time 7 < t; ; q@(tys) are nominal 


oF, 
I ierations and velocities at time f; ; and 5 (7) 
l 


is a Green’s function or ‘influence’? function which 
plays a decisive part in the theory of linearized trajec- 
tory perturbation mechanics [5], [7], [8], and |9]. 


Fixed Terminal Times and the Green’s Function 


In the next section of this discussion, Varying Termi- 
nal Time Defined by an External Constraint, the pur- 
poses served by the At, perturbation in Eq. (12) will 
become clear. However, in order to emphasize the sig- 
nificance of the Green’s function in linearized per- 
turbation theory, it will be assumed for the moment 
that the terminal time is fixed at ¢; for both nominal 
and off-nominal trajectories, i.e., Aty = 0. Then Eq. 
(12) reduces to: 


The Green’s function a (7), is therefore a second order 
l 


tensor operator which transforms the perturbation 
vector, Ag;(r), to a new vector whose components are 
parameter deviations, at the nominal terminal time 

4 The reader will note that in appearance and manner of 


derivation these equations are analogous to the differential 
correction equations of classical orbital mechanics [10]. 


a 


ty , from a standard trajectory. The partial differential 
coefficients which are the components of the influence 
tensor, may be termed, in the most general sense, 
fixed terminal time, first order perturbation coefficients. 
(Sometimes referred to as “error coefficients” in guid- 
ance and control studies.) Their evaluation will be- 
come a topic of investigation on the following pages. 
With reference to Eq. (13), expression (12) can now 
be written in another form, viz: 


Agu(ts) = Agu(ts) + G(tr) Aty . (14) 


Varying Terminal Time Defined by an Ex- 
ternal Constraint—The Modified Green’s 
Function 


The purpose served by the perturbation, At;, in 
the foregoing developments has perhaps been some- 
what obscure. It is clear, from the physics of the prob- 
lem, that At; , unlike the Agi(7)’s, does not actually alter 
a nominal trajectory. The only function of this “per- 
turbation” is to provide a means of comparing 
variables on the altered trajectory at time 


ty = ty + Aty 

to their respective nominal trajectory counterparts 
at time ty. When trajectory deviations are of interest, 
not at a fixed terminal time, but rather, at a variable 
terminal time determined by a certain criterion, a 
mathematical formulation which makes allowances 
for a Aty becomes eminently useful. The following ex- 
ample will attempt to demonstrate how Eq. (12) may 
be used to calculate Aq,(é;)’s, where the time ft; is 
now assigned to a particular event (in this instance, 
arrival at closest planetary approach) which, because 
of previously introduced disturbances, is displaced in 
time from ty; by a small increment At; . 

Consider a special case where the time ¢; corresponds 
to the point of closest planetary approach on some 
nominal interplanetary trajectory. Then it is reason- 
able to suspect that if Ag:(7) perturbations are intro- 
duced at time + < ¢;, the time ¢; will find the vehicle 
at some point other than closest planetary proximity 
on the off-nominal flight path. It follows, therefore, 
that if the terminal time ?¢; of interest is now the time 
of off-nominal closest. approach, then ?¢;, in general, 
differs from ty by a small increment At; which must be 
determined by means of the appropriate equation of 
constraint which must be satisfied at both the nominal 
and off-nominal times, ie., 7p(ts) = fp(iy) = 02 

The transformation equations which relate planet- 
centered spherical coordinates to inertial coordinates 
may be written in symbolic form which, for example, 
expresses the dependence of r, as 


rot) = Pla); vy), 2), a: (15) 


The time, ¢, appears separately here by virtue of the 
fact that the target planet would, in general, be moving 


5 Where r, is the vehicle-target position vector. 
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in some prescribed fashion against the background of 
inertial space. From Eq. (15), it follows that: 
dP dx , dP dy 4 OP oP dz , oP 


CEN 5e7 aan, Bea red (16) 


= Qlr(t), y(t), 2(t),4(t), y(t), 2() 


Since the parameters 2, y, 2, ete. are associated with 
the position of a vehicle, then Eq. (16) may be written 
in the previously established notation, viz: 


*p(ty) = Qlas(ts), ty] BS Me 6. 


Obviously, if the parameters q,(t;) and t; are subject 
to small variations, then 7,(t;) must also change. De- 
note by 7,(#;), the vehicle’s radial velocity on a per- 
turbed trajectory at time #; = ty + At, , and, by 7p(ty), 
the vehicle’s radial velocity on a nominal trajectory 
at time t;. Then 7,(#;) and 7,(t;) are related by: 


Fo(bp) = ty(ty) + Atel) = felts) + 
(17) 


I 8 (tat. 


$1 0 


Cs) Aga(ty) + — 


The significance of z (ts) in Eq. (17) is, in this in- 


stance, more forcefully demonstrated by a deductive 
reasoning process rather than straightforward manipu- 
lation of Eq. (16). 

Suppose At; = 0 in Eq. (17); then: 


Folly) = Tylty) + Ary te) = Fpl ty) 


+ EF Waele). (17a) 


Equation (17a) determines the off-nominal value of 
7, at the nominal terminal time t;. In general, an 
off-nominal radial acceleration, 7,(t;), will also exist 
at the same instant in time. Then, by linear approxi- 
mation: 


Ty (Ey) = Fp(ty) + Fp(tp)Aty = fp(ty) 


stn 2, - (ip)Age(ty) + Fy(ty)Aty. ye 


However, from Eq. (17a), it follows that: 


7 p(t) = (ty) =F is [Arp (ty) ] : 


d ; : f : 
But, a and A are commutative operations in this 
0 


instance. Therefore: 


fp(ty) = = 7, pkey) == Ar Grae 


Consequently, Eq. (18) reduces to: 


Ald) = Holt) + YS (taal) ae 
s=1 s a 


“1 [7 (ty) tr Af, (t;)] Al; . 
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Iti is clear that Fy(ts) is never zero, and AF (ty) Fisk in} 


initial disturbance is panne Therefore, is | 
first order of small quantities; 


6 fl 
7p(t,) = tp (ty) ae Le (t,)Agqs (ty) + y(t) Aty. (18 Db} 


Comparison of Eqs. (17) and (18b) leads to the f ol 
lowing conclusion: 


88 (4) = F(t) = OG)> 


OF « ; 
“dt 
Equation (18b) for its equivalent, Eq. (17)] m ys 
now be cast in still another form since an examinations 
of Eqs. (13) reveals that: - 


which parallels that arrived at for —— 


6 


OF, 
Aq.(t;) = De a 


l=1 


(r). (19)} 
Therefore Eq. (18b) goes over to: 


7p(ty) — ip(ty) ats Ss 2 ae (ty) 


s=1 l=1 Os 


OF, 
6qi 


7 )Aqi(r) + p(t, At, . 


The necessary conditions for r, to be a minimum on} 
the nominal trajectory at time f;, and on the off-| 
nominal trajectory at time fy , require that 7,(t;) = 0] 
and 7,(t;) = 0.° Under these circumstances Eq. (20)} 
reduces to: 


= 0+ My) ee 


t)Aty, COMI 


from which it follows that At; cannot be arbitrarily 
specified in the presence of a constraint, ‘as rather, IS) 
uniquely defined by the vector Agq:(7), i.e. 


=. Cr); Se | 
which constitutes an intuitively acceptable result, | 
since the displacement in time of closest approach is a 
direct consequence of the Aqi(7) perturbations. 
Admitting to Eq. (21a), the small changes Ag, (é;) | 
in the trajectory parameters at the point of off-nominal | 
closest planetary approach must be given by a special | 
form of Eq. (12), viz: 


At; = 


6 
= Vay 
Aq. (ty) = Se oe (7) Aqi(r) 
1=1 0q1 


(22) | 


a 6 6 74 
FOUTS) oat 9 
| 


* The vanishing of #, is a necessary but not sufficient con- | 
dition for rp to assume a minimum value. However, the dis- | 
tinction between maxima and minima can, in this example} 
be made quite readily through physical reasoning alone; | 


otherwise, the derivative #, , must be taken into considera. 
tion. 


Box (ty) = =e 3 | bn — vets) ou ) | 


1=1 s=1 Q(tr) 99s 


re Os mgs 
ie = (r)Aq(r); Osh = a a 


if Ants) = E = om i: w)|, then (22a) can 


be written in a more compact form given by: 


4 


all) = = = x A = (r)Aqi(r) 

E ; (22b) 
A = 2 An Aqi(r), 

F a 


where the second order tensor Aj: is simply a modified 
Green’ s function which is compatible with the mathe- 
‘matical constraint imposed at t; and 7, . 

_ One final word on this subject is in order: Although 
‘in the preceding example the mathematical constraint 
“was seen to enforce a stationary value for the vehicle- 
‘target position vector, the reader should bear in mind 
‘that a modified Green’s Function such as Aj, might 
easily be derived for a variety of terminal constraints. 
As a matter of record, the tensor A,(t;), as it appears 
in Eq. (22b), is representative of any pertinent mathe- 
matical constraint having the form Q(t;) = Q(i;) = 
0, (t; = ty + At;) and subject to the proviso that 
Q = Qla(t), d. 


Variable “Initial” Times 


It has previously been stated that the vector q;(t,) 
‘is a function of the vector gi(7) and the time t; counted 
from 7, 1.e.: 


gilts) = F(qi(r), ty). (23) 


However, if 7 is, in fact, an intermediate time and the 
set of six parameters gn(0) are associated with the 
time t = 0 S r St,, then: 


quT) = Gi(qm(0), ee (24) 
Consequently; 
au(ts) = FiiGilan(O), 7], ty}. (25) 


Equations (23), (24), and (25), together, are nothing 
more than a disguised form of the elementary integral 
relationship which guarantees that: 


ty OE ty 
fo atoa= [adnate [aly a. 06) 
0 0 T 
On the strength of Eq. (26), Eq. (25) can be carried 
one step further to: 
qu(ty) = Fr'[qm(O), ty] (27) 


and the intermediate time 7 completely disappears 
from the formulation provided that t; is now counted 
from t = 0. In this context, 7 retains its initial charac- 


ter only in the sense that a perturbation vector Aq:(r) 
is introduced at this point in time. 

An examination of the tensor A,, in Eq. (22b) 
discloses the dependence of its components on the 


nominal quantities q , 7», and a evaluated at the 


8 


nominal terminal time ¢;. But 7, and the derivatives 
aq) 
and q, , thus one is led to the rather remarkable result 
that the tensor A,, is uniquely defined by the standard 
trajectory’s initial conditions, its nominal terminal 
time, ¢;, and the nature of the terminal constraint, 
without consideration for the time of perturbation, 7+ or 
the perturbation vector itself. On the other hand, the 
magnitudes of the partial differential coefficients, which 


are themselves, readily derivable functions of 


are the components of the Green’s function me (7), 
1 


are related to the size of the contribution made by the 
second integral on the right-hand side of Eq. (26) 
since the Aqi(7) perturbation vector must exert its 
influence on q,(t;) through this integral. Consequently, 
as one might expect, these coefficients and the associated 
influence tensor are functions of the running variable 7; 
0 <7 St,. The situation is further aggravated by 
the fact that the functions /;, are, in general, not avail- 
able in closed form, therefore, the components of 
on 7) must be evaluated by a relatively laborious 
numerical procedure. However, in this instance, the 
end more than justifies the means, since a time history 
of the Green’s function (along with the invariant A 
tensor) completely describes the sensitivity of the 
q.(t;) variables to any arbitrary (small) perturbation 
vector introduced at any intermediate time 7 on a 
nominal trajectory. In fact, if care is taken to insure 
against a violation of the small disturbance assumption, 
the principle of superposition holds and the perturba- 
tion analysis can be extended to include situations 
where the components of Agi(7) are continuous (or 
piecewise continuous) functions of time on an interval 
T2 — 71, e@.g.: let SAgq(t;) denote the total increments 
in qx(ty); then: 


(t) 


(28) 


6 T2 
=D 4a] 2G) dale. 
al 

As an illustration of a situation where Eq. 28 is 
physically meaningful, one might reflect on a corrective 
guidance scheme which employs relatively low levels 
of thrust for finite periods of time (as opposed to im- 
pulsive corrections). Let w(t), 2 = 1, 2,3, represent 
the components of the imparted acceleration vector, 
T/m, where m is the mass of the spacecraft and T 
is the corrective thrust vector. Then; 


de (bys ="; ( 0 at; 
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and therefore, the total changes in terminal state af- 
fected by the control function are given by integrals of 
the form, 


mee CE say Hull dt, 


where [71 , 72] represents the interval during which the 
corrective thrust is applied. 
In any event, the evaluation of the Green’s Function 
_ (7) is, at once, the essence and central problem of 
i 
linearized perturbation mechanics. Attention is now 
directed toward two possible operational techniques 
which will generate the influence tensor. 


Evaluation of the Green’s Function by Direct 
and Indirect Integration 


From a conceptual standpoint, the simplest and most 
straightforward approach to be followed in evaluating 
the influence tensor is direct integration. Examine 
one of the thirty-six components of the Green’s func- 


tion, e.g., ae (7). By definition: 
1 


ony (7) = nee {Filqi(7) 
WANE Rone ead we 
aa Filqi(r), go(7), g3(7), t,]} /Aqi(r) 
or if Aq(7) is sufficiently small; 
5 (7) a (Pilar) 
+ Aq(r), g(r), g(r), +--+ , tl (29a) 
=e Filq(r), (7), w(t), tr]}/Aq(r). 


The second term in the numerator of Eq. (29a) is, 
with reference to Eq. (23) and the definitions of the q 
variables, the nominal x component of the vehicle’s 
velocity at the nominal terminal time ¢t; , and therefore 
is well known when the standard trajectory has been 
selected. The first term can be evaluated by means of 
direct simultaneous numerical integrations of the equa- 
tions of motion, which, for the restricted n-body prob- 
lem, take the form 


ie (hy ant), (30a) 

¥ = g(x, y, z, t), (30b) 
and 

2 = Ale,.y, 2 ty, (30c) 


where, once again, ¢ 
planetary movements. 

The integrations proceed from the time 7 at which 
the initial conditions (7) + Aq(r), Ag(r), ete. 
are impressed, to the nominal final time t; . Here g(r), 
q@2(7), etc. are nominal intermediate parameter values 


appears separately to allow for 
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available from the standard trajectory at time 7 ang 
Aq(r) is a small but otherwise arbitrarily selected 
variance. When the integrations have been accom) 
plished, the difference quotient in Eq. (29a) is com} 
pletely defined and the degree of accuracy of its ap) 
Dae to the Pa differential a 
OF 


Aqi(r) variance; the smaller the variance, the e } 
the precision. In a similar manner, the remaining 


thirty-five components of ele can be obtained. — 
LE: hon 


be eosin each time the lower limit of integration 
r is changed. Obviously, the difficulty is caused by thes 
varying values of qgi(r) and the changing interval off 
integration, ty; — 7, as 7 ranges from t = 0 to t = ty 
A complete time history of the Green’s function, 
generated by this technique, would therefore require a 
discouraging amount of numerical labor. F ortunately, 
another avenue of approach is available which achieves 
the same results as direct integration with a consider 
ably greater economy in computational steps. This: 
alternative method, which might be described as ani 
indirect integration, relies on the properties of thes 
adjoint to a system of simultaneous linear ordinary? 
differential equations and has been previously applied 
to problems of this kind by Pfeiffer [5], and to a related 
problem by Kelley, [7]. Supplementary discussions} 
may be found in Bliss, Goursat, and Goodman and] 
Lance, [9], [11], and [12]. Hoe in the interest of 
maintaining the self-sufficiency of this discussion, a) 
development of the necessary mathematics is included} 
in an appendix to this discussion. 

The trajectory Eqs. (30a), (30b), and (30c) take: 
the following form when linearized about a nominal | 
path: 


; 


esa erm of af of . 

AL =-— (Av), = Ae =. = H } 

ne (Az) ap 4. ay Ay + a Az, (Slag 

; 

Pe ee 0 ) 

age © (ag) = 0a a ee 

dt Ox oy dz (31b) | 

(At = 0}, | 

and | 
oh 

Aé (az) = Bas + Fay +o as, (31e) 


These three second order linear ordinary differential 
equations can be expanded to a set of six first order 
differential equations by employing the original defi- 


nition of the q variables, and the following identities: 
d ne eos ‘ 
i x) = Ad, (32a) 


d ' 
7 (Ay) = Ay, (32b) 


= (Az) = dé, (32¢) 


hen, Eqs. (31a), (31b), and (31c) can be written in a 
10rmal form, viz: 


d 6 
di (Agq.) = >. au Aq , (33) 
t=1 
where 
- Jo Te Re el | 
} Ox oy Zz 
3 7 
| Calarge Sob | 
> Ct aan eee, 
4 eae) 2 (at () oh oh oh | 
3 | dL) SOY +02 | 
4 | eth e (Oe 2 02.4,0 | 
Cee 08r 0 0: 70 
bs Cae Oa Oer 1) | 
. Meh ar of 
the derivatives aa etc., which are the time dependent 


‘coefficients of the linear equations, are evaluated along 
the nominal path. 

- The complete solution of Eq. (33) may be written 
‘in the following form: 


“ 


Aq(ts) = 2d, AniAgi(r). (34) 


At this point, it should be apparent that A,, above is no 
more than the influence function of the preceding dis- 
cussion. Furthermore, in accord with the results de- 
veloped in Appendix A, a complete time history of A,, , 
as 7 ranges from 0 to t; , can be generated by a back- 
wards integration of the system of differential equa- 
tions adjoint to Eq. (33). The adjoint equations are, 
in this instance; 


d : 
ae = a Dua AL; Dut = Ax - (35) 


d 


The solution vectors, \z(7), (ty; S 7 S 0) which cor- 


respond in turn to: 
1 0 
0 1 


0 0 
(ty) = | 0 Jace) = | 0 
0 0 
0 0 
are identically equal to the row vectors of the tensor 
A,,.' Thus the arduous process of direct integration, 


with sequential reinitializations is not really required: 
six uninterrupted backward integrations of the adjoint 


ete. 


7See Appendix A. 


equations with proper initial conditions imposed once 
and for all at ty achieve the identical result. 


Conclusion 


For the sake of brevity and generality, references to 
specific problems, which would necessitate explicit 
rather than functional equations, have been purposely 
avoided. In certain instances, the reader may have 
found this lack of specialized illustration an aid rather 
than a hindrance to understanding. The author does 
regret, however, that the protracted scope of this paper 
did not allow for a decidedly more complete discussion 
of magnitude limits for the perturbation vector com- 
ponents. The question of how small is a small per- 
turbation (within the framework of the first order 
assumptions) has not been adequately considered and, 
in all probability, cannot be completely answered out- 
side of a specific problem reference frame. Conse- 
quently, even though the function and operation of the 
linearized mechanics has been made reasonably clear 
in the main body of the text, the allowable extent of 
its application is still somewhat vague. 

Apart from all this, the amount of space devoted 
to the explanation of indirect integration deserves one 
last comment. When one appreciates the substantial 
amount of computing time required to derive a com- 
plete time history of the Green’s function by direct 
integration, the savings afforded by the indirect tech- 
nique are significant. Thus, from the point of view of 
economy alone, to say nothing of sophistication, in- 
direct integration is the more desirable approach. For 
this reason, the emphasis has been placed on the asso- 
ciated mathematics in order to provide a firm founda- 
tion for the application of this technique to specific 
problem areas. 


Appendix A 
Adjoint Differential Equations 


Consider the determinate set of n simultaneous, 
linear, ordinary differential equations in n dependent 
variables a; given by: 


oe = Gas ay ; a ilet (Al) 
dt eal 


where the a;;’s are real functions of the independent 
variable ¢ (which, for the purposes of this discussion, 
is time) and solutions are desired for the a,’s. 

The system of equations adjoint to Eq. (Al) may 
be written in the following form.” 


z (A) a => bi; Nae —— i n, (A2) 
dt j=l 


where 
big = Aji. (A2a) 
8 For the moment, the physical significance of the 2» 
variables is undefined. 
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It will now be shown that the integrals of the basic set of 
differential equations and the integrals of the adjoint 
set are interrelated in a very useful fashion. 

Since a; and \; are characteristically vectors in a 
parameter space, then their inner product is defined 
in the usual manner, viz: 


Sees (A3) 
i=l 


and the first derivative of this scalar function with 
respect to the independent variable ¢ is: 


d x d £ d 
=: . — i ; (h, ae a A4 
aon) Di Gat Lag (A4) 


which, with the aid of Eqs. (Al) and (A2) reduces to: 


= es a) = = Ss si Na Gh CH yy Ss (ony bi; dj ¢ (A5) 
jag= jh je 
But 7 and 7 are dummy indices of summation. They 
therefore may be interchanged in the second sum of 
Kq. (A5) to give: 


© ve) = Ee OO DED big het (CAG) 


i=1 j=1 t=1 j= 


But 6;; = a;; , therefore 
d 
— (rA-a) = 0. AZ 
7 (ha) = 0 (A7) 


Equation (A7) leads immediately to the interesting 
conclusion that the inner product of the a vector and 
the \ vector is invariant with time. Then, if the vectors 
a(t;) and \(t;) are associated with a specific value of 
the independent variable, it follows that: 

Dlr )a (7) = Do Ail ts alts), (A8) 
where 7 is any arbitrary value of ¢ different from t¢; . 

Suppose that in the adjoint system, the independent 
variable increases from t; in the direction of 7, iLe., 
ty <r. Then the A,(7)’s may be obtained by a simul- 
taneous integration of the Eqs. (A2) when a set of 
initial conditions \,(t;) for the adjoint variables has 
been specified. Let \;”” denote a particular set of ini- 
tial adjoint variables. Then Eq. (A8) adjusts itself to: 

DUM (a )aui(r ae (ty )a,(ty), (AQ) 
where the \{””(7) form the set of n adjoint variables at 
t = + which corresponds to an mth set of initial condi- 
tions df” (t+). In this context, Af”’(7r) and re” (ts) 
are still vector quantities. But if n cndependent but 
otherwise arbitrary sets of initial conditions, eG ee 
are considered as m ranges from 1 to n, then clearly, 
Ms” (r) and »§"(ty) can be completely contained by 
the second order tensors, \m;(7) and Nmilty). 

Then Eq. (A9) goes over to a set of n simultaneous 
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independent algebraic equations in a;(r) and a;(t,)} 
viz: 


S Nmi( 7 )ai(T) 


= > Ame(tyas(ty); m= 1, 1. 
i=1 ~ 


exact oe of that for the \ variables, ie., if 7 < ty | 
then the a;(7)’s are, in effect, initial” sanditee im) 
pressed at a varying lower limit of integration r 
However, the a;(ty)’s which correspond to these initia 
conditions may now be obtained by utilizing Eqs 
(A10) rather than by directly integrating the a Eqs: 
(Al). The advantage of this approach comes directly Lys 
from the fact that, because of the purposeful time flows 
reversal, the lower limit of integration, 7, for the 
system becomes the upper limit of integration for thes 
Ami(7) variables which, together with the assocratedt 
Ami(ty), and a;(r) are the coefficients of Eqs. (A10) 

Consequently, the \,,;(7)’s, and therefore the a;(t;)’s 
can be obtained by means of n simultaneous, unin 
terrupted backward integrations of the adjoint equa 
tions, where the mth integration is carried out for t es 
mth set of initial conditions contained in the Am:(t,)l 
tensor. ) 

In the most general cases (where \m:(t;) represents 
collection of independent, but otherwise arbitraril 
selected, sets of initial \ variables) Eqs. (A10) ares 
completely coupled, i.e., every a;(t;) appears in eachi 
equation. However, a judicious choice for the tensori 
Ami(ty) removes the difficulty, for if: 


on 


(0, mi 
Aras bp) = Oran = (All): 
1, m=7% 


then Eqs. (A10) simplify to: 


l,n (A12)) 
which are completely uncoupled algebraic expressions! 
for the an(ty)’s. Eq. (A12) is nothing more than the} 
readily recognizable standard form for the solution of} 
the set of linear homogeneous differential equations,, 
(Al). ) 


I 


DS Ami(7)ai(r) = amn(ty); m 
jal 
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A Study of the Martian Upper Atmosphere 


. and lonosphere’ 


Gilbert Yanow? 


Abstract 


- Hypothetical Martian models composed of 90-98% N» and 
10-2% CO» with varying surface temperatures and densities 
have been selected and a probable photochemical chain 
determined. This photochemistry has been solved for incre- 
mental layers from the altitudes of 650 to 250 km using the 
most recent solar data extrapolated to the mean distance of 
Mars. The results are presented in graphical form for density, 
temperature, mean molecular mass, electron density, and 
atmospheric gas and ion concentrations against altitude as a 
function of two and three body recombination. The results at 
present are preliminary, and while these data are admittedly 
speculation, it appears plausible to suggest the existence of a 
fairly dense, multi-layered, high altitude ionosphere about 
Mars, further suggesting that terrestrial communication 
techniques could be employed using, however, perhaps dif- 
ferent critical frequencies and skip distances. 


Introduction 


As the time when Man will actively explore the 
planets draws near, the need for upper atmosphere 
data of these bodies continually increases. Unfor- 
tunately, the amount of observed and experimental 
quantities of these data is limited. This situation has 
prompted the author to modify the astrophysical tech- 
niques used in stellar model calculations for determina- 
tion of atmospheric models. 

The method of approach is quite general in nature. A 
probable selection of possible chemical compositions 
for Mars has been taken from the findings of astro- 
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nomical observers. After a careful analysis, a chain of 
likely photochemical reactions for the selection was 
determined. Using the latest solar data available, the 
solutions to the photochemistry were integrated down 
through the media from 650 to 250 km. The parame- 
ters of surface conditions, recombination process, and 
chemical composition were varied to determine how 
critically they effect the photochemical solutions. 

It should be noted, that the numerical results are 
obviously quite speculative at the present time; they 
do, however, perhaps give an indication of which 
parameters are critical and some first order approxi- 
mation to them. 


Surface Conditions and Chemical Composition 


During the past years many observers have made 
estimates of the conditions on Mars. Recently de 
Vaucouleurs [1] has made an extensive survey of the 
field and determined the best values. This information 
appears in Tables 1 and 2. 

The surface pressure value can be derived by using 
the techniques of molecular light scattering. It is 
assumed that the atmospheric temperatures of the 
Earth and Mars are similar, and that the two at- 
mospheres have the same approximate scattering 
properties for a certain spectral region (A > 0.5 yw). 
The scattering coefficient ratio is then approximately 
equal to the gas mass ratio, and the surface pressure 
ratio is accordingly obtained by multiplying the mass 
ratio times the gravity ratio. The density value is 
directly obtained from these data. 
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TABLE 1 
Environmental Conditions of Mars 


surface pressure 85 + 4 mb = 4 X 10% dyne/cm? 
surface density = 10-4 gm/em? 

surface temperature = 273°K mean (sub-solar point) 
surface gravity 372.5 cm/sec? 


TABLE 2 
Martian Atmospheric Chemical Compositions 
Gas Volume (%) 
Ne 93.8 
Os 0.1 
A 4.0 (?) 
CO. =2.2 


The surface temperature may be calculated by con- 
sidering Mars as a gray sphere and applying Stefan’s 
law using the solar constant. Estimates of the surface 
temperature have also been obtained from radio- 
metric observations. The chemical composition, as 
given in Table 2, is based primarily on theoretical 
dynamics and chemistry. However, during observa- 
tions of 1947-48 Kuiper identified CO, bands in the 
near-infrared region of the planet’s spectrum. 


Solar Data 


Table 4 lists the majority of important solar lines 
in the region between 304 A and 1817 A. For the par- 
ticular cases of interest in this study, additional in- 
formation is given. The solution to the general problem 
is satisfied using the group of lines from 1808 through 
1400 A, the 1216 A line, the pair centered at 618 Ae 
the 584 A line, and the 304 A line. 


Photochemistry 


As a basis for calculation, a nitrogen-carbon dioxide 
atmosphere has been assumed. A possible chain of 
photochemistry for such an atmosphere is shown in 
Table 3. 

A rate coefficient of 10 *em*/sec [2] was used for 
reaction (1). The cross-section for this reaction was 
taken as 3.7 + 10°” em’ [3]. Values of the rate coeffi- 
cients for reactions (2a) and (6a) were not available. 
Consequently, the rather standard rate of 10 '’ em*/sec 
was adopted. For the reactions (2b), (5), and (6b) 
the temperature dependent relation 5 x 10“ 7" 
em*/sec [4] was assumed as a rate coefficient. C1t 
should be noted that the nitrogen atoms in reaction 
(6b) must be in an excited state; the corresponding 
reaction for atoms in the ground state is forbidden. ) 
The nitrogen molecules are assumed to follow a disso- 
ciative recombination process as is illustrated by reac- 
tions (3) and (4). The cross-section used for (3) was 
2X 10" em’ [5] and the rate coefficient used for 
(4) was 10° em*/sec [6]. The cross-section and rate 
coefficient applied in reaction (7) were 5 X 10° em? 
[7] and 10°’ cm*/sec [6] respectively. 
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TABLE 3 
Possible Photochemistry for a Nitrogen- ee 
Dioxide Atmosphere 


(1) CO: + hy— CO+0 
(2a) CO + O > COz + hy 
(2b) CO+ 0+ M-—CO.+ M 
(3) Ne+hy—Net+e 

(4) Nete—N’+N” 

6) N+0+M>NO7M 
(6a) N’ + N’—> No+ hy 
(6b) N+N+M—>N:+M 
(7) NO+hy->-N+0 i 
(8) CO. + hy ~ CO.* +e 
9) O+thyoOt+e 

(10) N+h-—N*+e 

(1) NO + hy ~ NOt +e 
(12) CO + hy + COt +e 


Reactions (8) through (12) describe the photochem- 
istry which determines the density and extent of th ; 
model Martian ionosphere. The cross-sections Be inh 


AT 


reactions (8) through (11) were, (8)3.3 x 10 


(Si 0). 12 10" em? [5]; (10) 10°-*’ em’; and a) 
2.2 x 10° cm’ [10]. The recombination coefficients 


employed for the reactions were, (8) 10’ cm’*/see} 
[2]; (9) 2 X 10°” em*/see [11]; and (11) 10 ° em*/see} 
TABLE 4 
Solar Emission Data [12] [13] 
Tine [Wave Lengtn|  Powet,Densigr | Borel Dau 
A erg/cm? sec quanta/cm®? sec 

He II 304 1 6.58 X 10° 
He I 584 0.05 6.32 X 108 
Mg X 610 0.05 7.55 X 108 
Mg X 625 0.01 1.51 X 108 
H Lyman « 938 0.004 — 
H Lyman 6 950 0.002-0.003 = 
H Lyman y 973 missing* a 
H Lyman 8 1026 0.06 Se 
Si III 1207 0.2 al 
H Lymana 1216 6 —— 
Si II 1265 0.02 —- 
Ol 1302 0.04 = 
Ol 1305 0.04 == 
Ont 1306 0.04 — 
Si II 1309 0.04 = 
Cor | 1335 0.3 — 
C Il 1336 0.3 == 
Si IV 1394 0.2 —= 
Si IV 1403 0.2 
Si II 1527 0.03 
Si II 1533 0.06 
CIV 1548 0.6 
(ON IA 1551 0.3 
CI 1560 0.12 5 & 10" 
He II 1640 0.3 
Ga 1656 0.3 
CI 1658 (On2 
Fe II 1670 0.2 
Si IT 1808 alc 
Si II 1817 ie 


* It is believed that H Lyman y is absorbed by the Earth’s 
atmosphere. 


[10]. The cross-section of reaction (12) was assumed 
to be the same as reaction (8), and the recombination 
coefficients of reaction (12) and (11) were assumed 
the same. 
_ Reaction (8) took precedence over reaction (1) 
when the wave length was equal to or less than 885 
. The same condition held for reaction (11) over 
f Bction ( ds) when the wave length equaled or was less 
than 1348 A. 
At present, values of rate coefficients, cross-sections, 
. recombination coefficients are not readily available. 
onsequently, some of the data used in the calcula- 
tions were not taken directly from the quoted reference, 
‘but may have been either an extrapolation or a mis 
Beeumed as characteristic of the particular reaction. 
No information was found for the spectral region of 
_approximately 500 A or less. It is known, however, 
that as the wave length approaches zero, the cross- 
“Section also goes to zero, and it appears Husihalle to 
assume a drop in the value of the cross-section of about 
one order-of- magnitude at the spectral region of the 
Hell line, 304 A. This assumption was made. 


Method of Solution 


* Assume that it were possible to insulate the model 
‘planet Mars from solar radiation. After a sufficient 
ength of time, all the electrons, ions, and atoms in 
the atmosphere would recombine to form only molec- 
ular nitrogen and carbon dioxide. Since for all prac- 
tical purposes the ionosphere photochemistry may be 
considered reversible, this equilibrium condition may 
be regarded as the state from which the ionosphere was 
formed. 


Initial Density Distribution 


It is assumed that the ‘insulated’? atmosphere be- 
haves isothermally, is in hydrostatic equilibrium, and 
the acceleration of gravity follows an inverse square 
relation; the density may be given by the customary 
density equation, 


po exp —(mgo/kT)(ah/a + h), C1) 


I 


p 
where 

po = surface density, 

m = mean molecular mass of the gas, 

go = surface acceleration of gravity, 

k = Boltzmann constant, 

T = abolute temperature, 

a = radius of the planet, and 


h altitude above the surface. 


It is further considered that in the upper regions of the 
model atmosphere the constituent gases will not be 
well mixed because of a low amount of convection, 
but will take on a distribution that is a function of their 


respective mean molecular mass, the gravity gradient, 
and the kinetic temperature. This distribution is ap- 
proximated by using Eq. (1) with the surface density, 
po, equal to the fraction of the constituent gas time the 
overall density. For example, if the gas were CO, 
and it amounted to 5% of the total atmosphere, the 
py used would equal 0.05 <X 10 * gm/cm’ for an over- 
all surface density of 10 * gm/cem’. 


Photochemical Analysis 


The number of photochemical reactions of any 
particular type may be expressed 


lim AiNnQr ) (2) 
where 
Qui = number of photochemical process occurring 
of specied 7 at altitude h, 
A; = cross-section of process of species 7, 


Nii = neutral particle density of species 7 at alti- 


tude h, and 
Q, = number of quanta of wavelength \ available 
per sec. 
Using Beer’s Law, for normal incidence 
AQ, = A:NniQn,Ah (3) 
and 
Q, = Qy, exp 2 > — ANpAh, (4) 


where Q,, = Q, at some datum level. 
A combination of Eqs. (4) and (1) yields 


qui = AN Qn exp 2, 2) — ANniAh. (5) 


Equation (5) can be applied directly to reactions (1) 
through (7). In reactions (8) through (12) electron- 
ion recombination is also a consideration. Bates and 
Massey [11] have shown that 


dN/dt = (q/1 + B) — (ae + Ba;)Ne 


(6) 
— (N./1 + d) dp/dé, 
where 
dN/dt = rate of change of the neutral particle, 
density or the rate of production of 
electrons, 
N, = electron density, and 
a,, a;, B = coefficients of various modes of recom- 
bination. 


Taking steady-state conditions, 1.e., 
dN/dt = d/dt = 
Kq. (6) reduces to 
q/1 + B = aNe, (7) 


where a = effective recombination coefficient. As a 
general rule 8 < 1, so Eq. (7) may be written 


Na= (Gfa)>: (8) 
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Equation (8) may be applied to reactions (8) through 
(12). The photochemical differential equation of 
reactions (1) and (2a) to determine the amount of 
CO and O at altitude h is given by 


a(CO)/dt = d(O)/dt (9) 
= Kyl, — KaCcoCo = 9, 
where 
K, , Koa = rate coefficients, 


vf 


ll 


photons or quanta absorbed per sec 
(I = qco,), and 
Coco, Co = concentrations of CO and O. 


The solution to this equation is 
Coo Ss Co = [Ki/Koa(N co,4,Qx)]”. (10) 


In a similar manner the solution for reactions (3), 
(4)> (5); and’ (6a) is 


Cy = { —K5Co€ x (11) 
PUA Kal C ead Qe fy 2a 


where C,, = third body concentration. The solutions 
to reactions (8) through (12) all have the form of 
Eq. (8). 


Boundary Conditions 


Mitra [14] has shown that an approximate value for 
the critical density or the density as a function of alti- 
tude where the atmospheric gases will just start to 
escape, 1s given by 


= [(mgoro/kT) — 2]/ro'r , (12) 
where 


molecular diameter, and 
ro = distance from the center of the planet to the 

point, where the atmosphere begins to be- 
have isothermally. 

Using the values of 

m = 28 X 1.67 X 10” gm 

go = 372.5 cm/sec’ 

ke = 1.38 610° erg /°K 

PES TVS007 1K 

ro = 3410 km (3400 km radius plus 10 km tropopause) 

5 X 10 ° em (after Mitra, 1952). 


o 


o 


The value of n., 4.7 X 10 "° gm/em’, obtained at the 
approximate altitude of 640 km, was adopted. Below 
the approximate region of 250 km recombination in- 
creases rapidly so that for all practical purposes the 
ionosphere has ceased to exist. 


Calculations 


The results presented in this paper have been cal- 
culated using a Bendix G-15 digital computer. The 
technique used is approximately the following. An 


initial density at a particular incremental layer, for 


example 645 km to 650 km, was found using Hq. (1). 
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Then using the energy of the 1216 A region, reaction 3 
(1) and (2a) were solved for the CO, O, and OG 
equilibrium condition concentrations. By successive? 
use of the 618 A, 584 A, and 304 A energies, reactions 
(30 (Ayan aid (6a) (or perhaps reaction, (3), (4),} 
(5), and (6b)) were solved for the N, N2, O, and NO 
concentrations. The 1808-1400 4 618 A, 584. A, 2 


ion-neutral equilibrium balance in reactions } 
through (12). When a photoreaction took place it was 


and in certain instances excite the atoms produced.) 
After these sequences occurred there was a certain} 
amount of energy in excess of what was used for each} 
event. For example, using the 618 A energy with reac-} 
tions (3) and (4) there was 0.7 eV extra per a y 
with 584 A 1.7 eV, and with 304 A 178 eV. This e 
cess energy must go into a kinetic form, and after ay 
sufficiently long period of time, this kinetic energy | 
should be fairly evenly distributed over the particles || 
of the particular incremental layer. A simple “3k7”7/ 
relationship was then used to calculated the increase 
in the temperature over the isothermal value produced 
by the photochemistry. The gas in the incremental 
layer was assumed to behave as a perfect gas and fur- 
ther that the temperature was the only parameter to 
change. Consequently, the effect that the temperature 
change had upon the density was to decrease the dens- 
ity by the ratio of the original temperature over new 
temperature. A new density was obtained in this man- 
ner, and the whole calculation of the photochemistry 
was repeated. 

Equation (4) was then used to determine the trans- 
mitted energy to the next lower layer, and using this 
energy the above procedure was carried through for | 
this lower increment. In this manner the problem was 
numerically integrated from 650 km to 250 km. 

: 


Results and Comments 


Figure (1) is a density vs altitude curve as a func- 
tion of recombination process. The curve is based on a 
surface density and temperature of 10° gm/em* and 
273°K respectively. From about 350 km to the lower 
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limit, the density value appears the same regardless 
of whether two body or three body recombination is 
used. However, in the upper regions the two body proc- 
ess produces a somewhat greater density. Also it is 


noticed that the two body curve does not display the 


discontinuity shown by the three body curve at the 
altitude of about 500 km. 

Figure (2) is similar to Fig. (1) with the only ex- 
ception being that the calculations were based on a 
surface temperature of 253°K rather than 273°K. 
The shapes of the curves of the two figures are the 
same. The only difference the twenty degree drop in 
the surface temperature makes is to lower the curve, 
point for point, about 50 km. When the surface tem- 
perature is raised, the curves are displaced upward by a 
similar amount. A change in the surface density has 
the same effect upon the solution as a change in the 
surface temperature. 

Figure (3) is a relation between temperature and 
altitude using the surface temperature and density of 
973°K and 10 * gm/cm’ respectively. Again the three 
body curve exhibits a rather large discontinuity at 
about 500 km, indicating temperatures that are per- 
haps unrealistic. The two body curve, on the other 
hand, represents values which are considerably more 
in the realm of reality. 

Figure (4) is also a temperature vs altitude graph 
but based on surface conditions of 253°K and 10° 


em/cem’. It shows the same characteristics as described 
in relation to Fig. (2). 

A curve of mean molecular mass against altitude is 
displayed in Fig. (5). From approximately 400 km 
downward both recombination methods appear to 
produce an atmosphere composed of primarily molec- 
ular nitrogen. Above this altitude, however, the com- 
position is strongly dependent on whether two body or 
three body recombination is employed. Using the lat- 
ter, the prime constituent is atomic nitrogen chile 
with the former the model atmosphere continues to: 
be composed of primarily molecular nitrogen. 

Figure (6) represents a profile of the model Martial 
ionosphere. Only data based on a chemical composi- 
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effect upon the solutions. The concentrations of N 
are almost identical to these curves, the only small fe TT 
departure occurring with the two body curve from 
about 400 km downward. As is demonstrated, the 
choice of recombination process strongly effects the ie | 
numerical results. t 
= 
The three body process produces an ionosphere of 8 
: ie ber 2 
two layers, separated by a marked discontinuity. E 
. s: < 
From an analogy with the Karth, the uppermost peak 200 ae ie ae ere ve 2 
was designated the 7, layer and the lower peak the F, Lu 
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sudes of about 300 km and 200 km with respective 
densities. of the order of 10° electron/em’ and 10° 
electron/em’. 
The two body model calculations also produce a 
two layer ionosphere. However, the layers are not 
nearly so obviously separate. These Ff, and F; peaks 
are at the elevations of about 430 km and 360 km, with 
density values of 3 X 10° electron/em’* and 4 X 10° 
-electron/ em® approximately. It should perhaps be 
noted that even when a recombination coefficient as 
large as 10 '° em*/see and a surface density of 10” 
gm / em® were used, peaks of the order of 10° electron/ 
ff em® and 10° electron/em® were still obtained for the 
calculated F, and F, regions. . 
a Figure (7) is a curve of atomic nitrogen concentra- 
¥ tion with altitude as a function of the recombination 
_ process. The results of this graph seem to indicate that 
3 in the regions of the atmosphere above 300 km, for the 
‘surface conditions used, two body recombination is far 
- more efficient than is three body recombination. 
Figures (8) and (9) are respective graphs of the 
atomic oxygen and carbon monoxide concentrations 
against altitude, depending on the recombination 
‘process employed. Of particular interest is the fact 
that the curves appear to indicate that the form of 
~ recombination is not of cardinal importance near the 
_ regions of the upper or lower limits of integration but is 
of importance in the mid regions of about 400 km 
altitude. 

Nitrogen oxide vs altitude as a function of recom- 
bination process is shown in Fig. (10). With a two 

_ body process the NO concentration is far less than with 
three body recombination. The two body NO layer 
is also much thinner than the three body layer, but 
both concentration curves peak at about the same 

altitude. 

Figure (11) is a curve of the nitrogen oxide ion con- 
centration. Only a three body recombination mode 
produces any reasonable amounts of NO*. 

The concentrations of the oxygen ion, depending on 
whether three body or two body recombination was 
used, are presented in Figs. (12) and (13) respectively. 
The two body process yields less O* than does the three 

body, and the peak of the two body curve occurs at 
an altitude of about 100 km lower than the correspond- 
ing three body curve. Further, while the two body 
process appeared to narrow the NO* concentration, 
the opposite is the case with the accounts of 0°. 


Conclusions 
Regardless of which recombination process was used, 
the two body or the three body, the results consistently 
indicated the possible existence of a multilayered, 
high altitude, fairly dense ionosphere about the planet 
Mars, if the atmosphere is composed of primarily 
nitrogen. The author feels that the results of the work 
show this technique of study promising, and is at pres- 
ent attempting to start the task of model Earth cal- 


culations. A better evaluation of the procedure can 
then be made with a comparison of the Earth study 
results with known terrestrial data. 
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Attitude Stability of an Elastic Body 


of Revolution in Space 


Leonard Meirovitch’ 


Abstract 


The motion of a body in space is described by three angular 
velocities called spin, precession, and nutation. 

A rigid body of revolution with no external moments has 
an angular momentum vector, h, constant in space and its 
spin or symmetry axis precesses about the vector h at a 
constant attitude angle @. An elastic body, however, dis- 
sipates energy and as a result the angle @ changes. If the 
decrease of angle @ denotes stability, it is shown on basis of 
energy considerations that stability is imphed by C/A > 1 
and instability by C/A < 1, where C and A are mass moments 
of inertia about the spin and pitch axes respectively. 

A body consisting of two elastic disks connected by a rigid 
shaft so that C/A < 1 is considered. Using the elastic solu- 
tion of a disk subjected to gyroscopic forces the time ¢ re- 
quired to reach an attitude angle 6 is obtained. 

The results are compared with other investigations of the 
same subject and the agreement is found very satisfactory. 


SYMBOLS 
h = angular momentum vector 
0 = attitude angle measured from vector h to spin axis 
A, C = moment of inertia about pitch and spin axes respec- 
tively 
U = strain energy 
Ye = hysteretic damping factor for normal stresses 
o = stress 
E = Young’s modulus 
V = volume 
fh = kinetic energy 
ey) = initial angular velocity 
t = time 
6 = nutation velocity 
D,;, = modified plate flexural rigidity 
r, a = polar coordinates 
p = mass density 
¢ = spin velocity 
D = plate flexural rigidity 
w = plate deflection 
y = Poisson’s ratio 
h = plate thickness 
W(r) = amplitude of plate deflection 
Yq = hysteretic damping factor for shearing stresses 
ty = cyclic time 
a,b = external and internal radii of the disk 
Introduction 


The motion of a moment-free rigid body with prin- 
cipal mass moments of inertia A, B, C is an unsteady 
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periodic precession and nutation about the resultant 
angular momentum vector h which is fixed in spac 
Steady motion is possible only about the axes of maxi 
mum or minimum moment of inertia. A rigid body of 
revolution A, A, C however has no a 
moment of inertia and, therefore, the moment-free 
motion consists of steady precession of the axis of sym- 
metry about the fixed angular momentum vector h. 
The angle @ between the symmetry or spin axis and the 
direction of the h vector is called attitude angle and is a 
constant [1]. 

A real body, however, is elastic. A vibrating elasti¢ 
body undergoes stresses and deformations resulting in 
energy dissipation which in turn causes a change im 
the attitude angle 6. The forces causing the elastic 
vibrations are the gyroscopic or inertia forces. The 
rate of energy dissipation per cycle of stress is the | 
area within the hysteresis loop as shown in Fig. 1. 

It is customary to assume that the energy dissi- 
pated in one cycle is a fraction of the strain energy. 
This is obtained by multiplying the strain energy ex- 
pression by a coefficient yz which is called the hys- 
teretic damping factor. 

In general, for a body subjected to stress of one kind 
only, namely co, the energy dissipation per cycle is 


given by 
vas = dn ffl em 
ee hee 5) a volume i se) 


The change in the attitude angle @ due to energy 
dissipation can be of a stabilizing nature or the op- 
posite. It follows that for an elastic body where energy 
is dissipated during vibration, the concept of stability 
of motion should be revised. It was shown [1], that the 
expression of the kinetic energy rate of change is 


dT end : dé 
eae Clan | eee S 108 6 — 9 
a Car @ 1) sin 6 cos 6 Wt (2) 
where 
A = mass moment of inertia about the pitch axis, 


; a : 3 : 
C’ = mass moment of inertia about the spin axis, and 
wy = initial angular velocity. 


Since energy is dissipated rather than gained, the 
kinetic energy must decrease. Hence, the right hand 


TARE SD Oe ee Oe Ge Ee ee 


Strain Energy 


Fig. 1. Stress-strain diagram showing the hysteresis loop 


) - precession 


Fig. 2. The body and the angular velocities 


side must be negative. There are two possible cases: 


G dé 
(a) gis ° and cee 
and (3) 
€ do 
(b) a 1<0O and 7 > 0; 


From Fig. 2, one observes that instability is charac- 
terized by a continuous increase of the initially small 
angle 6. 

Hence, 


C > A implies stability 
C < A implies instability. 


The above principle has been used in devices which 
stabilize satellites of the first type by intentional dis- 
sipation of energy. 

The problem to be discussed here treats the case 
shown in Fig. 2, where the two circular disks are elastic 
and the connecting shaft is assumed rigid. The elastic 
solution of the elastic disk subjected to gyroscopic 


forces was described in a previous paper [2]. Using the 
expression of the deflection, w, derived in the above 
paper, an expression for the strain energy is obtained. 
By equating the rate of energy dissipation, Uais , 
to the rate of change, 7’, of the kinetic energy, a rela- 
tion between the rate of change, 6, of the attitude angle, 
6, and the angle, 6, is obtained. Integrating the rate of 
change, 6, the time, ¢, to produce an attitude angle, 
6, is determined. 


The Kinetic and Strain Energy Expressions 


In the first reference, it was shown that the kinetic 
energy expression of a body of revolution is 


t= Scar |1 + ($= 1) sin? a]. 


The second reference shows that the deflection of the 
circular elastic disk subjected to gyroscopic forces is 
given in polar coordinates by 


(4) 


w= xa, LJiCyr) + X(2, 1) Yr) 


4X (SoD cry ER eee | © 


YDa 


-sin (a + ¢t). 


In the above Ji(yr) and Yi(yr) are Bessel functions 
of the first order and first and second kind respectively. 
Ii(yr) and K,(yr) are modified or hyperbolic Bessel 
functions of the first and second kind respectively. 
The coefficients X(j, 1) were determined from the 
boundary conditions. In addition 


y= ape (6) 
and 
f = pw e (2 — °) sin 6 cos 0. oe) 


The strain energy due to normal stresses is expressed 
in polar coordinates by 


2 
Ow 


= 5D ff | ay 20 — 
Us ae a4 Area | (vo) ae ”) or” 


Eh’ 
12(1 — »*) 
rigidity. Since every point at an angle a goes through 
all values ranging from —W(r) to W(r) twice in one 
revolution, it follows that the strain energy per cycle 
of stress is obtained by writing expression (8) in terms 
of the amplitude, W(r), of the deflection, w, and letting 
a vary between zero and 27. Consequently the strain 


where D = hDz = is the plate flexural 
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‘energy due to normal stresses per cycle of stress is 
Ustjevole = xD | {a'tx4, 1)Ji(yr) 
b 


4+ X(2,1)¥ilyr) — X83, hr) 
— X(4,1)Ki(yr)Pr — 201 — »)¥° 


[xan{-taon +| oun 


ian)) i X (2, 1) {-2 vr Yo(yr) 


2 f 
+|2:- | rin) + X(3, 1) {-+ Io(yr) 


9 (9) 
Ee Ee Bs | hian)} + X(4, 1) 


e Ki(yr) + Ee ue | Kw} 
| xa, 1) oan) a 2 rion} 
\ yr 


+x, {¥en — 2 rab} 

yr 

+X,0 {eon ~ 2 now} 
yr 


— X(4,1) {Kalan + 2 alan} |} dr 


In a similar way the expression for the strain energy 
due to shearing stresses per cycle of stress is written 


UO Gjoyete —e 2n(1 i v)D ipa 
{x4 1) | JaCar) cals Jin) | 
yr 
+x2,0[ Yr) - 2 nO) | (10) 
yr 
18 G@ehab) Ee Sigs hor) | 
yr 


— Xa 1) | Kula) + = Kw) |} dr dl 


Kinetic and Strain Energy Rates of Change 


The rate of change of the kinetic energy is obtained 
by simply differentiating Eq. (4) with respect to time 


yy ynomes C 9 A 
ID = (Cay (C - 1) sin 6 cos 66. eB) 


As mentioned in the first section, the energy dis- 
sipated in one cycle of stress is assumed to be a frac- 
tion of the strain energy. Specifically, 


BU jevcts == VnU zyeycle aii aU ejeycle 5 (12) 
where Usgyeyete and Ugyeycle are given by Eqs. (9) and 


(10), and yz and y« are hysteretic damping factors 
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for normal and shearing stresses respectively. To find 
the rate of change one has to divide expression (12) 
by the cyclic time f) which is the time required to per- | 
form one revolution. It was shown previously | 1] 

that the spin velocity expression is 


o= (1 = =) wo COS 6, 7 (139 | 
so that the cyclic time is 
rout or” 2a = ae 
ee soe: (1 2 (14 
= 4) antes 


It follows that the strain energy rate of change is 


i (1 = = wy COSA 
is i (ve Usjeyete + Ye U Stet (15) 


Attitude Angle as a Function of Time 


The rate of change of the attitude angle @ is obtainedl 
by equating the rate of change of the kinetic energy to 
the rate of energy dissipation. In order to account 
for the sign one writes 


a Ue, (16) 


and using expressions (11) and (15), Eq. (16) yields | 
rate of change of 6 


x il 
2a wo C sin 6 


(ye Onfazae Sy ae U ciemaene (17) 


6 - degrees 
Fic. 3. Rate of change 6 as function of attitude angle 6 


t — hours 


6 - degrees 


Fig. 4. Time ¢ as function of attitude angle 6 
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Note that Usyeycte and Ugjcyele are depending on 96. 
Equation (17) integrated gives the time ¢ as a func- 
tion of the attitude angle 6. Separating variables one 
obtains 


(18) 


Ulnievele = Ye Ugjeycie) 


where Usgyeycie and Uejeyele are given by Eqs. (9) and 
(10). 


= ‘ sin 6 
t = 2r wo G dé, 
0. (Yx 


Numerical Example 


Using the configuration shown in Fig. 2 a numerical 


_ example was worked out. Data used is as follows: 


Internal radius b = 4.0 in 
External radius a = 20,0in 
Thickness h = 0.5 in 
- Mass density p = 7.3237 X 1074 Ib 

A : in“ sec? 
Poisson’s ratio v = 0.30 
Young’s modulus KE =30X 10lbin-=? 
Hysteretic damping factor Yn, = 0.1 

for normal stresses 
Hysteretic damping factor vq = 0.1 

for shearing stresses 
Initial angular velocity w = 20 rad sec"! 

C/Ax—10st 


Equation (17) in conjunction with expressions (9) 
and (10) were used to determine the rate of change 
6 as a function of the attitude angle 6. The integrals 
in expressions (9) and (10) were evaluated numeri- 


- cally. The results are plotted in Fig. 3. In order to 


obtain the attitude angle, 6, as a function of the time, 
t, a procedure more suitable for computer use was 
adopted instead of equation (18). The function 6 vs. ¢ 
is exhibited in Fig. 4. 


Conclusions 


The results reported in the previous section indicate 
that the analysis presented is perfectly justified, and 
that the energy dissipation through cyclic stress is a 
factor which can become quite serious. 

In order to obtain the solution of the plate the cus- 
tomary assumptions of elastic behavior and small 
deformations were made. Both were justified by the 
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numerical results. In other cases, however, the assump- 
tions may not hold and large deflections and stresses 
beyond the elastic limit may occur. Though an analytic 
solution under such conditions does not seem possible, 
the fact remains that such a problem may exist. Hence, 
for cases C/A < 1 a careful analysis has to be made. 

The assumption that the nutation velocity, 6, is 
small compared to the spin velocity, ¢, and the pre- 
cession velocity, ¥, appeared to be fully justified. 

In the presently considered case, one can deduce 
from Fig. 4 that in a matter of hours the axis of sym- 
metry deviates from the initial direction along the 
angular momentum h to an appreciable angle 8. 

In general, the numerical results from the previous 
section conform quite well with the expectations. 
They also indicate excellent agreement with previous 
work done along the same line by Thomson as far as 
the dynamical behavior of the system is concerned. 
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Satellite Networks for Global 
Coverage 


Frank W. Gobetz’ 


Abstract 


This paper investigates the problem of establishing a net- 
work of satellites positioned such that at least one satellite is 
in view at all times from any point on Earth. In addition, the 
minimum number of satellites necessary to fulfill this condi- 
tion is determined over a wide range of satellite altitudes. 
Because of their symmetry two orbit patterns are considered 
for this application: distributions of polar orbits, and arrange- 
ments composed of orbits placed perpendicular to the faces 
of the regular polyhedra. The former are found to be superior 
in all cases. 

An analytical method is derived by which the optimum 
number of polar orbits for a given total number of satellites 
can be predicted. These optimized patterns are then con- 
sidered in detail and operating altitudes are determined for 
systems consisting of from 6 to 100 satellites. An important 
consideration in this phase of the analysis is the proper syn- 
chronization of satellites in adjacent orbits as well as within 
each orbit. Previous solutions to this problem have disre- 
garded such synchronization and are therefore unnecessarily 
redundant. 

Six satellites appear to be the absolute minimum number 
of satellites for the two orbit patterns considered. The opti- 
mum distribution of six satellites yields complete coverage at 
all times at an altitude of only 5260 n. mi. Some consideration 
has been given to the problem of obtaining fewer satellites 
with other patterns, and a special pattern requiring only five 
satellites is described. 

The calculated altitudes are found to compare favorably 
with ideal instantaneous coverage patterns composed of 
satellites placed at the vertices of regular polyhedra circum- 
scribed about the Earth. Although these instantaneous ar- 
rangements comprise the best physical distributions of satel- 
lites for obtaining global coverage, the relative positions of the 
satellites cannot be maintained over a complete revolution 
and the patterns are thus impractical for an operating system. 

Consideration is given to pattern stability and the polar 
array is shown to be an excellent choice from this standpoint. 
In addition, calculations are made to estimate the probable 
lifetime of each polar pattern, based on the satellite operating 
altitudes. 


Symbols 
ro Radius of Earth 
s Are-length radius of satellite coverage 
a Central angle of satellite coverage 
A Spherical area of satellite coverage 
h Altitude 
N Number of orbits 
n Number of satellites per orbit 
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n Total number of satellites 

y Angle formed by intersection of adjacen 
orbits 

Wo Spherical angle between orbits with co-diree- | 
tional motion 1 

y’ Spherical angle between orbits with opposing : 
motion z 


l “Tead;’” relationship between satellites in 
hie with those in an adjacent orbit 

Satellite positions measured from the North 
Pole 

m Satellite mass 

nN Lagrange Multipher 

Tr Any integer 

Z Arce length defined in Fig. 13 

L Satellite lifetime 

ro) Rate of rotation of the orbital plane about 


the Earth’s rotational axis 


D Orbital inclination to the equatorial plane 

le Period of revolution 

Jo Gravitational acceleration at the Earth’s 
surface 

S Satellite base area in drag calculation 

hs Satellite altitude based on restricted viewing 
angle 

Introduction 


}) 


Applications envisioned for artificial Earth satellites 


encompass such fields as meteorology, surveillance, 
navigation, and communication. A fundamental con- 
sideration inherent within each of these applications is 
the problem of distributing the satellites in such a way 
that every point on Earth is continuously in view of at 
least one satellite. The purpose of this investigation was 
to examine the problem of obtaining full-time, global 
coverage by a satellite network so conceived as to re- 
quire a minimum number of satellites. Although no 
attempt is made to prove rigorously that the distribu- 
tions discussed herein are optima, they nevertheless 
compare favorably with certain known ideal solutions. 


Assumptions 


All satellites are assumed to be at the same altitude 
and therefore only circular orbits are considered. The 
use of elliptic orbits apparently can not reduce the 
number of satellites required, and precession of the 


# perigee of such orbits might make them undesirable 
for this application. 
It was assumed that any number of satellites may be 
| uniformly distributed within a prescribed orbit and 
) that the orbits themselves may be attained precisely. 
Furthermore, all orbits in a given distribution contain 
_ the same number of satellites. This assumption is justi- 
- fied in part by Appendix I, in which arrangements with 
unequal numbers of satellites in each orbit are con- 
sidered. 
_ No stipulation is made as to which satellite will be 
visible to a given portion of the globe at a given instant. 
| One or more satellites will be visible at all times al- 
though the particular satellites in view will change 
_ periodically. This simplification is an important one in 
Be that it allows the Earth’s rotation to be neglected. 


~ Analysis 


PE 


The area on Earth visible to a satellite is a circle 
_ described by the locus of tangents from the satellite to 
| the Earth. For example, in Fig. 1 the satellite “sees” 
) the area between tangents SA and SB. Thus if s repre- 
sents the are radius of the small circle seen by the 
} satellite, 


Supe ek 
(e4 To 
8 ee (2) 
A s = = To 
therefore, a COs ae * a (3) 


It is apparent that the nondimensional radius, s/7ro , of 
this small circle is a function of the satellite’s altitude, 
_and that the circle approaches a great circle as a limit 
when the altitude is infinite. Thus the two quantities 
s/ry and h/ro can be used interchangeably for a given 
Earth radius 7). For convenience, s/79 will be used in 
most of the diagrams and in much of the explanatory 
discussion which follows. It is important to realize 
that the radius s/ro represents the altitude of a satellite 
in any diagram, even though the satellite itself may not 

_ be shown. 

Since the assumption has been made that each orbit 
will contain the same number of satellites, it can be 
shown that the total number of satellites in any pattern 
may not be less than six (see Appendix IT). That is, no 
orbital pattern can maintain global coverage for the 
entire period of revolution if less than six satellites are 
involved. A total of six could be achieved either by two 
satellites in each of three orbits or by three satellites in 
each of two orbits. Similarly a total of eight would entail 
two in four or four in two, etc. 

With regard to the placing of the orbits, two arrange- 
ments are considered here: (1) an arbitrary number of 
orbits which intersect at two common points 180 deg 
apart; due to the Earth’s oblateness perturbation the 
intersection points would normally have to be at the 


Fic. 1. Coverage of one Satellite 


poles, but this will subsequently be discussed in more 
detail, and (2) orbits oriented so that they will be 
parallel to the faces of the regular polyhedra. For sim- 
plicity these two patterns will be referred to as “‘polar- 
type” and “‘polyhedron-type”’ distributions. 

In determining altitude requirements for the various 
satellite patterns to be considered here it will be neces- 
sary to establish “critical” satellite orientations as 
design conditions. Obviously when the projections of 
two satellites on the Earth are in close proximity, the 
region between them is well covered. But if three satel- 
lite projections are regarded as the vertices of a spheri- 
cal triangle ABC as in Fig. 3a, coverage of the midpoint 
0 of the triangle requires that each satellite be at the 
altitude corresponding to s/r>o. The problem then, is 
one of finding the most serious “hole” opened up by the 
relative motions of the satellites, and establishing a 
minimum altitude requirement for its closure or cover- 


age. 


Polar Orbits 


In analyzing the polar-type orbit arrangements the 
critical satellite orientations are easily ascertained from 
Fig. 2. Since the satellites are uniformly spaced in each 
orbit the most critical separation of satellites occurs 
where the orbits are most distant from each other. This 
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Fia. 2. Polar orbit paths 


condition occurs 90 deg from the orbit intersection 
points, i.e., at the equator. Conversely, the most redun- 
dant coverage occurs over the intersection points where 
satellites move closer together. 

Assuming that satellites in two adjacent orbits move 
in the same direction, the critical condition for covering 
the region between the two orbits is prescribed by some 
position of the spherical triangle ABC in Fig. 3a. Im- 
agine a circle centered at 0 and passing through points 
A, B, and C. The radius, s/r¢, of this circle is maximized 
by moving the three vertices of the triangle in their 
respective orbits, always maintaining synchronization. 
This maximum value of s/r> represents the required 
altitude for coverage between two orbits in which the 
satellites move in the same direction. Unfortunately 
such synchronization is not possible between all ad- 
jacent orbits in a given polar-type pattern. The mo- 
tions of satellites in at least two adjacent orbits will 
always have opposing sense and the altitude require- 
ments thus imposed must be at least as bad as those 
found with co-directional motion. 

The critical satellite orientation when the satellites in 
adjacent orbits move in opposite directions is one which 
exhibits symmetry about the centerline between orbits 
as in Fig. 3b. The determination of satellite altitude 
necessary to cover the region within the four points D, 
I, F, and G entails circumscribing a circle about the 
spherical quadrilateral formed by these points. The 
radius of the circumscribed circle is some s/ry which 
represents a higher altitude than that found when the 
motion was synchronized. Thus a condition such as that 
shown in Fig. 3b determines the operating altitude of a 
system of satellites in a polar-type orbit arrangement. 

As an extension of this line of analysis, and for reasons 
which will subsequently be made clear, it is convenient 
to define an upper bound to the necessary satellite 
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S/fo= RADIUS OF COVERAGE OF ONE SATELLITE 


c) 


Fig. 3. Altitude requirements for polar orbits 


altitude. This upper limit can be established by disre- 
garding the positions of individual satellites and con- 
sidering only a ‘‘band-width” covered by the satellites 
in any orbit. A simple determination of this bandwidth 
is illustrated in Fig. 3c. Two satellites J and K in the 
same orbit are in positions which are symmetric with 
respect to the “equator.” If the magnitude of s/ry is 
such that both satellites can “see” point P, which lies 
on the equator and midway between the two orbits in 
the diagram, then the bandwidth extends one half the 
distance to the next adjacent orbit. Since each orbit 
contains the same number of equally-spaced satellites, 
the bandwidths are equal and the condition defined in 
Fig. 3c is more than sufficient to insure global coverage. 
That is, the magnitude of s/ro represents an altitude, 
Nmax » Which is always equal to or greater than the 
previously defined operating altitude. 


: 
| 
: 
: 


_ This upper altitude limit can be expressed as a simple 
unction of the number of orbits, N, and the number of 
satellites per orbit, 7. Similarly the altitude required to 
_ cover the region between two orbits whose satellites are 
| synchronized is a simple function of N and n. Since the 
latter altitude is always less than the operating altitude 
/ and the former is always greater, a definite bound is 
| placed on the desired operating altitude. 

' Determination of the operating altitude involves find- 
| ing the midpoint of a spherical triangle whose sides are 
_ known. While this is not a particularly difficult opera- 
tion, it involves the solution of several equations con- 
taining natural and inverse trigonometric functions 
which are too unwieldy to manipulate in an optimiza- 
tion operation. The use of the upper and lower bounds 
will thus be useful in the optimization procedure of 
- Appendix III. From this analysis it becomes possible to 
predict, for a given total number of satellites, , the 
“number of orbits and satellites per orbit which will 
yield a minimum altitude solution. In this way the num- 
ber of cases which must be examined in detail can be 
_ reduced, and at the same time limits can be established 
' within which the desired altitude must fall. 
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As yet no stipulation has been made as to the rela- 
tionship between satellites in adjacent orbits. Since the 
spacing of satellites within each orbit is uniform, satel- 
lites in adjacent orbits with co-directional motion will 
maintain the same phase relationship or ‘lead’ to 
each other. A logical choice for good synchronization 
might be a lead of 3 the spacing between satellites in a 
given orbit. Another distribution which was investi- 
gated is that in which the lead is 27/n?. A procedure by 
which the best lead may be rapidly determined is out- 
lined in Appendix IV. 


Orbits Parallel to the Faces of Polyhedra 


Configurations possessing a high degree of symmetry 
can be created if the satellite orbits are constructed 
parallel to the faces of the regular polyhedra. The orbits 
of each such network intersect at regular intervals thus 
describing intricate patterns of equilateral spherical 
polygons as shown in Fig. 4. 

Placing a satellite at each intersection, the directions 
of motion can be adjusted so that collisions will never 
occur. Because the polygons formed by the orbits have 
equal sides, it can be concluded that the same sym- 


() 


CUBE 


ICOSAHE DRON 


% IDENTICAL ARRANGEMENTS OCCUR 
FOR THESE TWO POLYHEDRA 


Fic. 4. Polyhedron orbits 
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metrical pattern will be repeated after every Qr/n 
radians of revolution. Furthermore the condition of 
maximum satellite separation must occur when each 
satellite is midway between two interesection points, 
and consequently it is this position which determines 
the proper altitude for global coverage. 


Discussion of Satellite Distribution 
Polar-Type Distributions 


When the polar orbit distributions are optimized 
according to Appendix III, the two curves of Fig. 5 can 
be constructed. The optimum polar configurations are 
then easily determined for any number of satellites. 
For example, when » = 6 one could choose either 
N = 2orN = 3, but since the former lies closer to the 
curves it is the better choice. Similarly, when y = 32, 
N = 4, and so on for all values of 7. 

Certain arrays are not valuable however, even though 
they represent the best solutions for given 7’s. One such 
case is » = 16, for which Fig. 5 predicts NV = 2 rather 


than N = 4. But when the minimum satellite altitude 
“ WT 
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Fig. 5. Optimum polar configurations 
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Fia. 6. Polar view of orbit spacing 


is calculated, it is found to be greater than that for 
n = 15, N = 3. There are several such redundant cases, 
of which a few are 7 = 16, 22, 25, 27, 30, --- etc. | 

The small circles in Fig. 5 renee optimized polar 
distributions which are not redundant. That is, each 
increase in the number of satellites entails a reduction, — 
however slight, of the required altitude. 

It has been previously stated that for polar arrays, 
the limiting satellite orientation is prescribed by oppos- 
ing motion in adjacent orbits, and that only two such 
orbits need exist regardless of the total number of 
orbits involved. Obviously the altitude requirements can 
be reduced if the orbits are spaced so as to favor this 
particular pair, the technique being illustrated in Fig. 6. 

When the orbit spacing is optimized in this manner, 
considerable reductions in altitude are achieved. Table 
I is a summary of the optimum polar distributions, in- 
cluding the appropriate leads l, optimized angles Wout . 
and the minimum altitudes, Aopt . 


Polyhedron-Type Distributions 


Table II summarizes the altitude requirements for the 
polyhedron arrangements. The column headed hmin 
refers to the symmetrical positions with satellites located 
at the intersections of the orbit paths, while the operat- 
ing altitudes based on the worst satellite orientations 
are shown in the last column. None of the polyhedron 


TABLE I 
Characteristics of Optimized Polar Arrays 
' 
4 I 
F : : * (deg) ts es 
6 2 3 40.0* 77.1 5260 
8 4 45.0 83.3 2860 
10 5 14.37* 86.7 2340 
12 3 4 ODD 46.3 1798 
15 5 72.0 49.0 1194 
18 6 10.0* 55.4 1037 
21 7 2050 53.6 840 
24 + 6 30.0 31.2 732 
28 i 7.34* 40.1 618 
32 8 22.5 37.3 486 
40 5 8 5.62* 31.4 425 
45 9 20.0 28.9 326 
55 11 16.35 30.5 275 
60 6 10 18.0 217 241 
72 12 15.0 24.3 201 
84 7 12 Baliye Up \h 194 
91 13 13.83 19.6 156 
105 15 12.0 21.3 137 
112 8 14 12.85 16.7 124 


* Denotes | = 2z/n?; all other leads are | = r/N. 


TABLE II 
Polyhedron Distributions for Full-Time Coverage 
Figure N n n ania aia 
Cube 3 2 6 2515 oo 
Tetrahedron 4 3 12 1425 3645 
Dodecahedron 6 5 30 610 743 
Icosahedron 10 9 90 565 609 
TABLE III 
Altitude Summary for Instantaneous Polyhedron Distributions 
Figure n h/ro Description 
Tetrahedron Anil e2 40 Four equilateral triangles 
Cube 8 | 0.732 | Six squares 
Octahedron 6 | 0.732 | Hight equilateral triangles 
Dodecahedron | 20 | 0.237 Twelve regular pentagons 
Icosahedron 12 | 0.320 | Twenty equilateral triangles 


The regular star polyhedra are excluded from this summary 
because only the convex polyhedra represent minimum altitude 
solutions. Star polyhedra have 12 vertices and all are inferior 
to the regular icosahedron from the standpoint of altitude. 
(see Refs. 6 and fH) 


distributions was found to be better than a polar ar- 
rangement requiring the same total number of satellites. 


Comparison of Distributions 

In order to properly evaluate the merit of the solu- 
tions arrived at here, a comparison was made with the 
absolute minimum altitudes as defined in Appendix V. 
This appendix deals with a fictitious case in which the 
satellite coverage areas are “distorted” from their 
circular shape in order to reduce the overlap to zero. 
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Fic. 7. Instantaneous polyhedron arrangements 


While these solutions are not physically attainable their 
usefulness lies in the fact that they represent a condition 
of zero redundancy. 

Another arrangement useful for purposes of compari- 
son is one composed of satellites placed at the vertices 
of a regular polyhedron circumscribed about a spherical 
Earth. Due to its symmetry this configuration is the 
best possible physical solution but provides only pe- 
riodic coverage. That is, the orbits cannot be so placed 
as to retain the positions of the satellites relative to each 
other and thereby maintain continuous coverage. A 
summary of the altitude requirements for each of the 
five regular polyhedra is provided in Table III and Fig. 
7 depicts the satellite arrays. 

Figure 8 compares these two idealized cases with the 
real solutions of Tables I and I and illustrates several 
interesting points. Both the smooth curve and the 
lower step curve approach an infinite number of satel- 
lites when the altitude is close to zero. However as the 
altitude is increased the ideal curve smoothly ap- 
proaches two satellites while the real solution remains at 
six until infinity is reached, at which time the final step 
in the curve goes to two satellites. As expected, the 
solutions corresponding to the instantaneous regular 
polyhedra fall in between the two curves. Due to the 
conclusions of Appendix I there is no advantage to be 
gained by filling the orbits with unequal numbers of 
satellites if 7 > 5, and therefore the lower step curve of 
Fig. 8 represents the best results obtained in this 
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analysis. As discussed in Appendix I, the case of five 
satellites in three orbits probably does not constitute a 
better solution, since the operating altitude is several 
times that required for six satellites. 

The satellite operating altitudes discussed thus far 
have been based on the assumption that a satellite can 
“see” an area on Earth described by the locus of tan- 
gents from the satellite to the Earth as in Fig. 1. Even 
if the Earth were a smooth sphere this assumption would 
be an optimistic one since the Earth’s atmosphere 
causes a refraction of light and radio waves. This effect 
is most pronounced near the “twilight”? region of the 
satellite coverage area because an observer in this 
region must look through a greater thickness of at- 
mosphere in viewing the satellite. 

To compensate for this effect the right angles SAO 
and SBO in Fig. 1 must be increased by some angle, 6. 
As a result, the satellite altitude must be increased to a 
value greater than h in Fig. 1. The resulting altitude, 
hs, is plotted against the ideal operating altitude A in 
Fig. 9. Therefore, if the angle 6 is assumed to be zero, 
h and hs are equivalent, but for 6 > 0, hs > h. 


Network Stability 


The solutions labeled “optimum” share the common 
feature that they depend upon a complex synchroniza- 
tion of satellites and orbits. Such precision is to no avail 
however, if this synchronization cannot be maintained 


120 The Journal of the Astronautical Sciences 


for reasons of stability. At least a qualitative considera- 
tion of perturbations to which the satellites will be sub- 
jected must therefore be undertaken. | 
Satellite perturbations arise from three main sources: 
the Earth’s oblateness, the atmosphere, and thir¢ 
bodies such as the Sun and Moon. Of these, the first 
two are by far the most significant (Refs. 3 and 4), | 
and accordingly third-body effects will be disregarded | 
here. Solar radiation pressure could also cause a sub- | 
stantial perturbation if the satellite surface area to mass 
ratio were unusually high, as in the case of the Kcho } 
balloons. Since this study is not directly concerned wit 
passive communication satellites however, this per- | 
turbation will be neglected. 
The oblateness of the Earth has a threefold effect o 
satellite orbits. 
1. Rotation of the major axis within the orbital plane | 
2. Reduction of the period : 
3. Rotation of the orbital plane about the Earth’s | 
axis c 
Since only circular orbits have been considered here 
rotation of the major axis is not pertinent. Effects on 
period of revolution cannot be disregarded, but these 
may be compensated by adjustment of the altitude of | 
the satellite orbit. 
Rotation of the orbital plane is predicted by the ex- 
pression (see Ref. 2) 


3.5 
= ro . deg 
@ = 10 — cos 7 i (4) 


It can be seen that the rate of rotation is directly related 
to orbital inclination, 7. Thus, polar orbits are unaf- | 
fected in this respect, and circular orbits of equal in- | 
clination are all affected identically. This is quite im- | 
portant in the present application, since the positions of 
satellites relative to each other, are not relative to the 
Earth, are of primary concern. 

With this in mind, a careful re-examination of the 
two-orbit case (V = 2) reveals an important exception. 
Regardless of the location of the orbit intersection 
points, the two orbits can have equal inclinations to the 
equatorial plane. Therefore if the intersection points 
are not at the poles, the orbits rotate at the same rate 
about the Earth’s axis and satellite coverage remains un- 
changed. This is a characteristic of only the two-orbit 
solutions, however, and is not true for more than two 
orbits. 

Satellite perturbations due to the Earth’s atmosphere 
result in a decrease of satellite altitude which advances 
rapidly as the dense region near the Earth is approached. 
Thus a pattern consisting of 78 satellites in a very low 
orbit has a lifetime of only a few weeks, while a six- 
satellite pattern at several thousand miles altitude can 
be expected to remain in orbit indefinitely. Satellite 
lifetimes are plotted against 7 in Fig. 10. In the prepara- 
tion of this graph the method of Ref. 2 was followed 
and numerical data was obtained from Refs. 4 and 5. 

It can be concluded with reasonable assurance from 
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TABLE IV 
Unequal Numbers of Satellites in Each Orbit 
N a b c d n 
2 2 4 = = 6 
3 4 == = if 
3 5 = == 8 
3 6 aa a 9 
4 5 “= = 9 
4 6 oe = 10 
5 6 = = 11 
5 a = — 12 
3 1 2 2 — 5 
2 2 3 == 7 
3 3 4 — 10 
2 3 3 = 8 
3 4 4 == 11 
4 2 2 3 53 10 
3 3 4 4 14 


the preceding stability considerations that polar orbit 
patterns are the best practical arrays for global cover- 
age. It had been previously stated that symmetry 
alone was a good argument for the use of polar distribu- 
tions, but the inherent stability of these orbits consti- 
tutes another important feature not likely to be found 
in other arrays. 


Concluding Remarks 


An important consideration in this analysis was the 
proper synchronization of satellites in adjacent orbits as 
well as within each orbit. Previous attempts at solution 
of the global coverage problem (Refs. 1 and 8), have 
disregarded such synchronization and are therefore un- 
necessarily redundant. 

Since only two orbital patterns have been investi- 
gated, the results presented here cannot be considered as 
absolute optima. However, a fundamental premise of 
this analysis has been that the desired optima will 
possess a high degree of symmetry. Figure 8 shows 
clearly that the instantaneous polyhedron distributions, 
which are extremely symmetrical, fall close to the step 
curve, and one can therefore conclude that if improve- 
ments are to be made they must at best be small. 
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Appendix I 
Unequal Numbers of Satellites in each Orbit 


When each orbit is required to contain the same 
number of satellites, a restriction is placed on the tota 
number of satellites in the pattern. Thus if 


= nN (5 


where 


distributions containing » = 5, 7, 11, 13, 17, 19, ete. | 
are excluded. In order to obtain such values of 7 it is” 
necessary to relax this restriction. | 

The cases investigated in detail are summarized in | 
Table IV in which a, b, c, d, --- represent the number | 
of satellites in each orbit wherea Sb Sc Sd---ete. — 
Each distribution was examined only until it could be 
shown to be inferior to a corresponding case in Table I 
or Table II and therefore the altitude requirements are 
not listed. 

It should be pointed out that these results are incon- | 
clusive in one respect. When adjacent orbits contain un- 
equal numbers of satellites and the satellites in each 
orbit are spaced uniformly, there is a considerable lati- 
tude of choice in selecting an appropriate lead between 
satellites in adjacent orbits. Due to the trial-and-error 
nature of this phase of the analysis, only a few possibili- 
ties could be examined in each case, and the results are 
therefore uncertain. In most cases however, the alti- 
tudes calculated were prohibitively great and could 
not be made competitive by any juggling of the satellite 
leads. 

Only one instance was found in which it was advan- 
tageous to fill the orbits with unequal numbers of satel- 
lites. This is the five-satellite distribution which is 
composed of three orbits. The method of coverage in 
this case takes advantage of the basic asymmetry in 
any pattern composed of an odd number of satellites. 
In Fig. 1la two satellites are shown in the same orbit in 
such relative positions that they leave a lune-shaped 
area uncovered. This uncovered area, only half of which 
is visible in the diagram, rotates in a clockwise manner 
as the satellites revolve in the plane of the paper. Be- 
cause two finite-altitude satellites cannot cover a great 
circle (see Appendix II), the center of the hemisphere in 
Fig. lla is within the shaded area. Two satellites can 
also be placed in an orbit which is perpendicular to the 
first as shown in Fig. 11b. Again only half of the region 
left uncovered by these two satellites is visible in the 
diagram, and the sense of rotation of the shaded area 
proceeds from the lower part of the diagram toward the 
top. 
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Fig. 11. Global coverage by five satellites 


If these two cases are how combined as in Fig. Ile, 
the entire sphere is covered for the instantaneous satel- 
lite positions shown in the sketch. When the satellites 
are allowed to revolve in their assigned directions, how- 
ever, the shaded regions soon meet, indicating a loss in 
coverage. Thus in Fig. 11d coverage 1s first interrupted 
at point D and at another point D’ which is in the hemi- 
sphere hidden in Fig. 11d. Note that the arc length 
D-D’ is only slightly greater than 90 deg. The un- 
covered (cross-hatched) region then spreads until in 
Fig. lle its area reaches a maximum. Then, as in Fig. 
11f it decreases in size to two points E and FE’ which are 
analogous to D and D’ in Fig. 11d. By visualizing the 
further rotation of the satellites it is not difficult to see 
that the uncovered regions arise at symmetrical posi- 
tions on the globe (180 deg apart). Therefore a single 
satellite could be placed in a third orbit and syn- 
chronized in such a way that it could cover these regions 
as they arise. 

Due to the complexity of this system a general oper- 
ating altitude has not been determined, although the 
altitude must exceed that for an instantaneous tetra- 
hedron distribution, which is two Earth radii. One 
special case was examined in detail and the result is 
probably close to the minimum altitude which can be 


obtained with a five satellite network. In this system the 
plane of the fifth satellite’s orbit would be viewed on 
edge in Fig. 11 and the satellite would pass from left to 
right on the visible side. For this particular case the 
minimum operating altitude was found to be 13 Earth 
radii or 44,720 n. mi. 


Appendix IT 
Minimum Number of Satellites 


In the body of this report the assumption is made 
that each orbit will contain the same number of satel- 
lites, n. As a consequence of this assumption distribu- 
tions containing less than a total of six satellites cannot 
achieve full-time global coverage. The following is a 
proof of this hypothesis. 

1. Consider two satellites at infinite altitude as shown 
in Fig. 12a. Each satellite is in view of a great circle 
(i.e., a hemisphere) and, together the two satellites 
maintain continuous global coverage. 

2. Tf these two satellites are now brought in slightly 
from infinity to some very great altitude a condition 
such as shown in Fig. 12b results. A band of width e re- 
mains uncovered because forh < © a satellite cannot 
cover a great circle. In fact this statement can be ex- 
tended to include any two satellites. That is, two satel- 
lites at finite altitudes cannot cover a great circle. 

3 Consider a large number of finite-altitude satel- 
lites in a single orbit as in Fig. 12c. Two circles of di- 
ameter ¢ remain uncovered. It can therefore be con- 
cluded that for h < «, at least two orbits are required 
for global coverage. 

4. Because the two uncovered regions in Fig. 12¢ 
are 180 degrees apart, they lie on a great circle. In step 
2 it was shown that at least three satellites are required 
to cover a great circle. Therefore, if N = 2,n 2 3, and 
n = 6.Since N > 1 by step 3, it can be concluded that 
n 6. 


IV IV 


Appendix 1 
Polar Pattern Optimization 


When the total number of satellites in a polar orbit 
pattern is specified, the exact distribution of satellites 
and orbits is not uniquely determined. For example, an 
arrangement consisting of twelve satellites may be com- 
posed of two, three, four, or six orbits with six, four, 
three, or two satellites in each. When the number of 
satellites, n, becomes large the numerous possible distri- 
butions become quite prohibitive to any trial-and-error 
approach. Clearly, a technique is desired by which this 
choice can be made less laborious. 

Due to the complex synchronization of satellites in 
adjacent orbits, the desired altitude requirements are 
not reducible to a simple analytical form. However, the 
altitude limits Amax and Amin MAY be expressed as rela- 
tively simple functions of N and n. If these limits are 
minimized holding » constant, the optimum distribu- 
tions should be defined by the solutions thus obtained. 
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(a) 


(b) 


Fig. 12. Proof of minimum number of satellites 


Referring to Fig. 3c, the maximum altitude limit can 
be expressed by: 


T 
COS — COS - 


To ia} 2N 
If the total number of satellites is held fixed, 
n(n, N) = Constant (7) 


Rmax can be minimized as follows: 
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Solutions of Eqs. (6) through (9) results in 


Tv Tv 
tan 4 - tan oN aa 
n 2N 
therefore, n = 2N, (11) 


Liiders has also obtained this result in Ref. 8. 

A similar procedure leads to optimization of hmin . It 
is necessary in this case to make the assumption of a 4 
lead between satellites in adjacent orbits, but this re- 
striction is not severe. Thus satellite C in Fig. 3a is 
assumed to lie on the equator and 


Rewin me 


—1. (12) 


It is expedient here to rearrange Eq. (12) to a different 
form. Thus 


wT Tv 
a BOC cae Oe 
F(n, N) = tan = = —* Lionel’ 


0 


sin — 
N 
The parameter F(n, N) may be minimized in place of 
Rinin since it is the are radius s/ro which is really being 
optimized in either case. As before 


oF 0 
d he Ase 
pun ant aN (15) 


Solution of Eqs. (13) through (15) with the constraint 
of Eq. (7) yields 


5 Tv Tv Tv T 
sin ta = en Qc 
N n COs = COS =. (16) 


Although this expression cannot be solved explicitly 
for either variable, a graphical solution is easily ob- 
tained. This graphical solution together with the solu- 
tion for Amax in Eq. (11) is shown in Fig. 5 as curves of 
_vs N. In order to determine N for some particular y, it 

is necessary only to choose the N closest to the curves. 


Appendix IV 
Optimum Lead 


__ The lead is of considerable importance because it 

affects not only Amin, but also the critical satellite 
orientation illustrated in Fig. 3b. The least desirable 
condition is one which possesses symmetry about the 
equator as well as about the centerline between orbits. 
A judicious choice of the lead can prevent this limiting 
orientation, thus directly lowering the required altitude 
for complete coverage. 

Consider the pattern pictured in Fig. 13a where one 
complete orbit and half of each remaining orbit are 
shown. For convenience the satellite placed at the 
upper intersection belongs to the completed orbit. 

The desired limiting condition involves the shaded 
area since the orbits bounding this region have opposing 
motion. If the lead is l, a satellite in the interior orbit 
must be at 


y = UN — J), (17) 
where y is measured from the upper intersection point. 


Similarly a satellite in the exterior orbit is at 


y= HO. (18) 
n 


If these satellites are now allowed to move until they 


are in the adjacent positions depicted in Fig. 1b, 


a) 


b) 


Fic. 13. Determining the optimum lead 


Also 
Qn 
Zryt+ aes (20) 
But such a condition will occur after each 2/n radians 


of revolution. Thus 
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where 

Five 0) eyo .ce (22) 
Given n, N, and 1 it remains now to determine the value 
of I. for which the satellites assume positions most 


nearly symmetric with respect to the equator. This con- 
dition is represented by a minimal value of the quantity 


\Z—y| =7(n-4—2r— 7 (N -1)) (23) 


nr 
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This, of course, is subject to the condition that T take 
on only integer values. 

The worst possible lead in any case would be one for 
which |Z — y|min = 0 and conversely the optimum 
lead is one for which | Z — y |min iS 2 Maximum. 

Equation (23) is perfectly general in that it facilitates 
the examination of any lead and the determination of a 
characteristic value of the quantity |Z — y|min. To 
compare two leads, | Z — y |min is calculated for each. 
The larger value corresponds to the better choice of 
leads, in the sense that it will yield a lower satellite 
operating altitude. 

A wide choice of leads was not considered in arriving 
at the results shown in Table I. It was found that leads 
of r/n and 27/n? were close to optimum in several cases 
and in all subsequent calculations only these two leads 
were considered. The relatively small reductions in 
satellite altitude which may be obtained by investiga- 
tion of other leads does not appear to justify the exten- 
sive calculations which would be required. 


Appendix V 
Perfect Distribution of Area 


As shown in Fig. 1 the area visible to a satellite is cir- 
cular in shape. Therefore if the Earth is to be covered 
entirely, considerable overlap of adjacent circles must} 
take place, resulting in redundant coverage. The area of) 
the zone cut out on the sphere by a single satellite is de- 
termined as follows: (Referring to Fig. 1) 


a T5— hy ra fo (24) . 
OE ene | 
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If the areas covered by each satellite were perfectly 
distributed with no overlap, the number of satellites re~ 


quired for global coverage would be 4 


z 
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Book Reviews 


Review of Methods of Celestial Mechanics, by Dirk Brouwer 
and G. M. Clemence. 

At the present time the interest in celestial mechanics, 
particularly in its application to space vehicles, is in the midst 
of a remarkable revival. In spite of this rebirth of interest, 
the new textbooks available have been few in number; and 
those that were available were weak—weak particularly with 
respect to emphasis on numerical applications. It is, therefore, 
with considerable delight that students of celestial mechanics 
welcome the new book by Dirk Brouwer and Gerald M. 
Clemence, entitled Methods of Celestial Mechanics (Academic 
Press, 1961, $10.50). 

Although Methods of Celestial Mechanics does not really 
consider applications to space vehicles, i.e., astrodynamics 
(only Section 14 of Chapter XVII specifically treats Earth 
satellites), and concerns itself very little with orbit determina- 
tion, it does lead the reader clearly and logically down the 
road to a thoroughgoing understanding of orbit prediction 
techniques. Both special and general perturbations are 
treated in considerable detail, as are the subjects of aberra- 
tion and the method of least squares. 

In Chapter I are derived the equations of relative motion 
and the vis-viva equation, as well as a number of other rela- 
tionships useful in the study of elliptical motion. Section 21 
introduces matrix notation and calls attention to the original 
work of Cayley, while the exceedingly useful f and g series 
are developed in Section 30. 

Fourier series and harmonic analysis are presented in 
Chapter II, and a useful definition of general perturbations 
is included; ie., “. . .expansion in Fourier series where the 
arguments are retained as literal quantities while the coeffi- 
cients are expressed as numbers.” In order to lay down a 
fundamental background for the analysis of orbits about an 
aspherical planet, Chapter III develops the gravitational 
attraction between bodies of finite dimensions. In particular 
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pages 125 ff. exhibit formulas that are well suited for the 
analysis of the motion of an object about a triaxial ellipsoid 
(such as the Moon). 

Chapter IV discusses the algebraic, tabular, and graphical 
methods for the representation of data and then proceeds 
with an excellent presentation of polynomial interpolation 
formulas. Lagrange’s (including a generalization), Everett’s, 
Bessel’s, and Newton’s (forward and backward difference) 
formulas are all discussed and well illustrated with numerical 
examples. The way is thereby logically prepared for numerical 
integration. Here they show, among other things, that the 
“yrobable error of numerical integration after n steps is 
0.1124n3/2 (in units of the last decimal), ... .”’ Certain useful 
symbolie operators are also defined in this chapter. 

In Chapter V Cowell’s method is applied to the close ap- 
proach of Ceres on January 6, 1941 and a derivation and 
numerical example of Encke’s method are given as well. 
With respect to the choice of special vs. general perturbations 
the following words bear considerable significance: “It cannot 
be known yet whether either method will gain at the expense 
of the other (numerical vs. analytical), but it is certain that 
the practical celestial mechanician will always profit by the 
use of a judicious combination of numerical and analytic 
methods.” 

The following chapter includes a lucid discussion and 
definition of planetary and stellar aberration. The next 
chapter (Chapter VII) then discusses the comparisons of 
observation with theory, dwells briefly on a combination of 
special and general perturbations, and concludes with good 
discussions of precession, nutation, and geocentric parallax. 
Of some note here is the equation for sidereal time referred 
to the equinox of 1950.0 found on page 208. 

A thoroughgoing discussion of error theory, including ac- 
cidental and systematic errors, Gaussian distribution, singu- 
larities, and ways and means to solve by least squares are 


iscussed in Chapter VIII. Chapter IX comes as close as 
ny to the subject of orbit determination in its treatment of 
ifferential correction, although the reader is really given no 
nsight into the question of the generation of a preliminary 
bit. This problem is presently one of the central ones in 
strodynamics, if not modern celestial mechanics. 

General integrals and equilibrium solutions are considered 
n Chapter X, where the force function, F, is introduced. 
The Jacobian integral, Tisserand’s criterion, surfaces of zero 
relative velocity, and the Lagrangian points are all discussed. 
In this latter regard they utilize the standard method for the 
treatment of small oscillations in mechanical problems; i.e., 
they develop linear differential equations with constant 
coefficients by expanding the coordinates of an object in the 
neighborhood of a libration point in a Taylor’s series and 
retain only the first-order terms. For a Trojan minor planet 
they find that the motion about the equilibrium position 
consists of a libration with a period of about 147 years and 
a more rapid superimposed oscillation with a period slightly 
in excess of the period of Jupiter (about 11 years). 

The chapter on variation of arbitrary constants is a rather 
involved one and requires careful study. Lagrange’s brackets, 
Delaunay’s method, etc., are all discussed. The lunar theory, 
which is usually understood to be “the problem of finding 
the Moon’s motion under the gravitational attraction of the 


arth and Sun, all three bodies being treated as mass points,” 


is thoroughly investigated in Chapter XII. It was of particu- 
lar interest to the present author that Brouwer and Clemence 
note that Brown’s tables of the Moon don’t reflect the best 
values of the lunar parallax and mass. The evection (largest 
periodic perturbation in the Moon’s longitude—about 1° 16’), 


“the annual equation, and the parallactic inequality are all 


discussed. It is also of interest to repeat the comment that 
none of the other classical methods can compete with Rabe’s 
dynamical method for evaluating the solar parallax. Hill’s 
Junar method and Brown’s differential correction procedure 
then conclude the chapter. 

The perturbations of the coordinates, which seems rather 


analogous to Encke’s method and effectively and interestingly 
utilizes the vis-viva equation, is the subject of Chapter XIII. 


This same chapter introduces the reader to Hansen’s device 
and Brouwer’s special perturbations method. Hansen’s 


“method (Chapter XIV) is felt to be the best procedure for 


treating first-order perturbations because its rapidly con- 
verging series 1s applicable to higher orders of 7 and e than 
most other methods. Furthermore, “Hansen appreciated the 
advantages of applying all perturbations of both long and 
short period to the . . . mean anomaly.” An interesting treat- 
ment of mean and osculating elements can be found in Section 
8 of this chapter. 

Chapter XVI demonstrates how to obtain secular terms 


-by utilizing Lagrange’s method. Section 6 specifically dis- 


cusses Jacobi’s determinant solution in which one reduces 
the off-diagonal terms to zero (a procedure that has been 
much used in differential correction of late). Other topics 
such as the limits on é and tan ¢ in the solar system and 
families of minor planets are introduced here also. The con- 
eluding chapter on canonical variables investigates infinitesi- 
mal contact transformations and other aspects of 
Hamiltonian mechanics. Here we find the only application 
to artificial satellites. Of particular interest in this regard is 
an exposition of Hori’s procedure, in which he removes the 
singularity for i near 63° 26’ by expanding in terms of the 
square root of the second harmonic. 

From cover to cover Methods of Celestial Mechanics is an 
excellent book. It is not truly a book that one without some 
experience or prior introduction to celestial mechanics or 
astrodynamics can easily fathom—but it does evolve each 
topic from first principles and should be comprehendible to 


an. advanced graduate student from other scientific disci- 
plines. Certainly for the graduate celestial mechanic or 
astrodynamicist it is a book well worth studying and repre- 
sents a notable contribution to the field. 
Rosert M. L. Baker, JR. 
Department of Engineering, 
University of California, 
Los Angeles, and 
Lockheed California Company 
Astrodynamic Research Center 


Review of The Physical Principles of Astronautics, by 
Arthur I. Berman. 

The Physical Principles of Astronautics by Arthur I. 
Berman (John Wiley and Sons, N. Y., $10.50) is a book 
designed to give the college student a concise and thorough 
exposition of the basic principles of astronautics, “. . .the 
science and technology of the placement and control of 
manned and unmanned objects in space.”’ The book is based 
upon a course instituted by Professor Berman in 1958 and 
develops most of the topics covered from first principles. 

Those possessing a background in modern physics, such 
as Berman himself, will find the book to be of particular 
advantage. Chapter One discusses the principles of astronomy 
and divides the field into cosmogony, astrophysics, and 
celestial mechanics. Here as in the rest of the text the areas 
of practical astronomy, navigational astronomy, and the 
burgeoning field of meteoritics have been neglected. These 
disciplines enrich the background of astronautics and should 
probably have been treated, if only in a cursory fashion. 
Aside from these omissions, the introductory material is clear 
and presented in an exceedingly interesting fashion. 

The physics of the solar system chapter is equally well 
done and a careful reading shows that the author has made 
every effort to include the most up-to-date data and astro- 
dynamic constants available at the time of publication—a 
very painstaking task indeed. If any criticism is to be leveled 
at the subsequent chapter on the foundations of mechanics, 
it would be that it is a bit too clear, too elementary. It would, 
of course, serve as a valuable review but hopefully the reader 
would have had some prior exposure to such subjects as 
Newton’s laws. On the other hand the discussion of accelerated 
reference coordinates is very relevant and deserves the elab- 
oration given this chapter. In dealing with potential and 
kinetic energy the vis-viva equation is developed and utilized 
in some simple illustrative examples. It is felt, however, that 
the full power of this exceedingly useful little equation is 
mitigated by the numerous different ways in which it is 
formulated. 

The subject of the dynamics of flight is also effectively 
presented but suffers somewhat from the overly emphasized 
analogy with potential barriers so useful to atomic physics 
and the lack of discussions of orbit determination and orbit 
prediction techniques. Another drawback can be found in 
the otherwise laudable attempt to utilize exclusively the 
metric system of units. Accuracy, simplicity ef formula, 
and historical precedent all recommend the dimensionless set 
of astrodynamic units (e.g., the Earth radius, a.u., speed at unit 
distance, etc.). The chapter on dynamics of flight, as well as 
many of the other chapters, is well illustrated by numerical 


such examples are indispensable tools in class 


examples 
instruction. 

Transfer orbits and space navigation are clearly developed, 
again with a number of useful numerical examples. The only 
two criticisms that apply here are that “A fundamental rule 
of astronautics” (p. 179) is not clearly shown to be a logical 
consequence of the variation of parameters, and that the use 
of differential correction is not discussed. The chapter on 
orbit perturbation 1s well illustrated and gives the student 
a clear intuition into the physies of “what's going on.” 
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Finally, a chapter on propulsion dynamics concludes the 
book. This chapter is particularly well done and shows that 
the author is knowledgeable in almost all areas of conventional 
and unconventional propulsion. Our only regret is that so 
interesting and well done a topic was not extended even 
further 
From almost every point of view The Physical Principles 
of Astronautics is a fine book and should be recommended 
as an introductory text to most fields of astronautics. In 
fact as Professor Berman states in the preface, ‘“This book . . . 


will serve to lay the groundwork at an intermediate level, fo 


works in astrodynamics . . . ,” and presumably also for works 


in propulsion, spacecraft structures, aerodynamics, communt- 


cations, and the other divisions of astronautics. 


Rosert M. L. Baker, JR.) 
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_ The Editors will appreciate the cooperation of 
uthors in using the following directions for the prepara- 
jon of manuscripts. These directions have been com- 
diled with a view toward eliminating unnecessary cor- 
‘spondence, avoiding the return of papers for changes, 
ind reducing the charges made for “author’s correc- 
ions.” 

Manuscripts 

4 Papers should be submitted in original typewriting 
(if possible) on one side only of white paper sheets, and 
should be double or triple spaced with wide margins. 
However, good quality reproduced copies (e.g. multi- 
lith) are acceptable, An additional copy of the paper 
will facilitate review. 


Company Reports 

_ The paper should not be merely a company report. 
If such a report is to be used as the basis for the paper, 
appropriate changes should be made in the title page. 
Lists of figures, tables of contents, and distribution 
lists should all be deleted. 


Titles 


The title should be brief, but express adequately the 
subject of the paper. A footnote reference to the title 


should indicate any meeting at which the paper has 


been presented. The name and initials of the author 


should be written as he prefers; all titles and degrees or 


honors will be omitted. The name of the organization 


with which the author is associated should be given in a 
separate line to follow his name. 


Abstracts 


An abstract should be provided, preceding the intro- 


duction, covering contents of the paper. It should not 
exceed 200 words. 


Headings 


The paper can be divided into principal sections as 
appropriate. Headings or paragraphs are not numbered. 


| Tllustrations 


Drawings should be made with black India ink on 
white paper or tracing cloth, and should be at least 
double the desired size of the cut. Each figure number 
should be marked with soft pencil in the margin or on 
the back of the drawing. The width of the lines of such 
drawings and the size of the lettering must allow for the 
necessary reduction. Reproducible glossy photographs 
are acceptable. However, drawings which are unsuitable 
for reproduction will be returned to the author for re- 


- drawing. Legends accompanying the drawings should 


be typewritten on a separate sheet, properly identified. 


Security Clearance 

Authors are responsible for the security clearance by 
an appropriate agency of the material contained in the 
papers. 


Mathematical Work 

As far as possible, formulas should be typewritten. 
Greek letters and other symbols not available on the 
typewriter should be carefully inserted in ink. Each 
such symbol should be identified unambiguously the 
first time it appears. The distinction between capital 
and lower-case letters should be clearly shown. Avoid 
confusion between zero (0) and the letter 0; between 
the numeral (1), the letter |, and the prime (’); between 
alpha and a, kappa and k, mu and w, nu and », eta and n. 

The level of subscripts and exponents, should be 
clearly indicated. Vectors will be set in bold face type. 
Authors should indicate this in their manuscripts by 4 
wavy underscore. 


Greek Alphabet 


A a alpha (a) N v nu (n) 

B £6 beta (b) ci xi (x) 

T y gamma (g) O o  omicron (8) 

A 6 delta (d) Il wr pi (p) 

E ¢« epsilon (@) P p_ rho (r) 

Z ¢ zeta (z) = os sigma  (s) 

H 7» eta (6) T 7 tau (t) 

6 ¢@ theta (th) > 2 upsilon (u) 

I «iota (i) ® ¢¢g phi (ph) 

K «x kappa  (k) X x. chi (ch) 

A »~ lambda (1) vw y psi (ps) 

Myp mu (m) Q w omega (6) 
The Orbital Elements 

a = semimajor axis 

e = eccentricity 

Q = longitude of the ascending node 

i = inclination 

«® = argument of perifocus 

T = time of perifocal passage 


Complicated exponents and subscripts should be 
avoided when possible to represent by a special symbol. 

‘Fractions in the body of the text and fractions occur- 
ring in the numerators or denominators of fractions 
should be written with the solidus as follows: 


cos (7x/2b) . 
cos (a/2b) () 


The intended grouping of handwritten formulas can 
be made clear by slight variations in spacing, but this 
procedure is not acceptable in printed formulas. To 
avoid misunderstanding, the order of symbols should 
therefore be carefully considered. Thus: 


(a + bx) cos t cos t (a + bx). 


In handwritten formulas the size of the braces, 
brackets, and parentheses can vary more widely than 
in print. Particular attention should be paid to the 
proper use of braces, brackets, and parentheses 
(which should be used in this order). Thus: 


{[a + (b + x)] cos ky}? (2) 


is required rather than ((a + (b + cx)") cos ky)?. 
Equations are centered, numbered, punctuated as 
complete sentences; and referred to in text as (2). 


is preferable to 
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References should be grouped together in a bibliog- 
raphy at the end of the manuscript. References to the 
bibliography should be made by numerals between 
square brackets [4]. 

The following examples show the approved arrange- 


ments: 

for books—{1] BAKER, R.M. L., Jr.and MAKEMSON, 
M. W., An Introduction to Astrodynamics, Academic 
Press, New York, 1st ed., 1960. 

for periodicals—{2] LAMORE, LEWIS, “Celestial 
Observations for Space Navigation,”’ J. Astronaut. Sci., 
6 (1959), 1-10. 
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